Mass added due to gravitational potential in a simple galaxy model.

gravity_mass.py

import random

# Gravitational constant
G = 6.674e-11 # m³ / (kg s²)

# Gravitational constant times the mas of the sun (nominal solar mass parameter)
GM = 1.3271244e20 # m³ / s²

# Sun's radius
r = 6.957e8 # m

# Sun's mass
M = GM/G # kg

# Speed of light
c = 299792458.0 # m / s

# Milky way diameter
mw_length = 9460730472580800.0 * 150000 # m

# Added mass to the system when separating the stars by the given distance:
def pair_mass(dist):
	# Potential grativitational energy when separating two sun-like
	# stars by the given distance
	pair_energy = (GM*M) * (1.0/(3.0/4.0 * 2**(1.0/3.0) * r) - 1.0/dist)
	return pair_energy / c**2

# Sample a point in the galaxy disk:
def rand_disk_point():
	while True:
		x = random.uniform(-1.0, 1.0)
		y = random.uniform(-1.0, 1.0)

		if (x*x + y*y) <= 1.0:
			return (x * mw_length * 0.5, y * mw_length * 0.5)

# Sample a star distance in the galaxy disk:
def rand_dist():
	a = rand_disk_point()
	b = rand_disk_point()
	return ((a[0] - b[0])**2 + (a[1] - b[1])**2)**0.5

star_count = 9e10

# Visible mass of the galaxy
lm_mass = star_count * M

print(lm_mass)

count = 0
accum_pair_mass = 0.0

# I assumed that sampling with different star distances
# and averaging (as a Monte Carlo simulation) would improve
# the result over the time, but it turns out the final distance
# is irrelevant.
#
# The initial distance is much more important. To improve the simulation
# to accurately describe milky way, I would need to sample stars' mass and radius
# from a accurate distribution of star types in the milky way, and would need
# a good model to tell what radius they would have if merged,
# (and black holes would be excluded, for obvious reasons).
while True:
	count += 1
	# Distance has almost zero effect of the pair mass...
	# It may be necessary to sample star masses from real
	# star distribution.
	pm = pair_mass(1.2848827903755582e+21) #rand_dist())
	accum_pair_mass += pm

	# Average potential gravitational energy mass accumulated
	dm_mass = (accum_pair_mass / count) * star_count * (star_count - 1)

	print("Added mass:", dm_mass, "Ratio:", dm_mass / (dm_mass + lm_mass))

	# Breaking because there is no point in carrying a Monte Carlo
	# simulation where the result doesn't change.
	break