Generating Point clouds using the Random Sampling function on the Navarro Frenk White profile for further simulation of the forces acting inbetween the stars. The end goal ist to simulate a stable rotating galaxy.

This is the [en] version of this post, the [ger] version can be read here

The goal of my Project is to generate realistic galaxies and dark-matter-halos.
In order to do this, I'm using the so called "Navarro-Frenk-White" (NFW) Profile that in combination with the Random-Sampling method can be used to generate the position of stars in a coordinate system.
By comparing simulated galaxies with real galaxies, a difference can be noticed: the stars get pulled away from a non visible object. This has been explained using the theory of "dark-matter", matter that is not visible. Its position and size can only be estimated by looking at its effects e.g. the path of the stars in the sky compared to simulations.

The NFW profile is a profile used for simulating the mass distribution in N-Body simulations. In a simplified way, the profile gets the distance $r$ from the middle of the galaxy and returns the probability that the star is at the given position. The profile is built up in the following way:

$$\rho_{NFW}(r) = \frac{ 1 }{ \sqrt{ 2 \pi } \cdot \sigma } \cdot \exp \left( \frac{ -\phi(r) }{ \sigma^{ 2 } } \right)$$

$$\phi(r) = \frac{ 4\pi \cdot G \cdot f_{0} \cdot R_{s}^3 }{ r } \cdot ln{ \left( 1 + \frac{ r }{ R_{s} } \right) }$$

Ok, don't panic, it's not as complicated as it looks. The complete function is actually not that interesting. The main focus lies on using it, so this is all we actually need to care about:

$$\rho(r) = \ldots$$

So we input $r$ and get a value that we can use. Let's do this:
We've got a star with the coordinates $x_1$, $y_1$ and $z_1$ and want to find out what the probability is that the star is generated. We can input the coordinates into the 3D Pythagorean Theorem:

$$r = \sqrt{x_1^2 + y_1^2 + z_1^2}$$

The result can be put into the NFW-profile resulting in a value $s$:

$$\rho(r) = \ldots = s$$

Now we've got a value $s$ and we want to find out if the star is going to be kept or if another star should be generated. We can to this by calculating $a = \rho(r_{min})$ and $b = \rho(r_{max})$. After generating a random value $r$ in the interval $[a;b]$, we can test if our value $s$ is greater or smaller than $r$:

If $r > s$, the stat should not be generated (it just isn't there) and if $r < s$, the star is generated (it's there).

After doing this for a lot of stars (for about 45000 original stars, on gets kept and the others are just not there), the following image can be seen:

Pseudocode for generating n stars:

# Generate a single star
def gen_star():
star_found = False

while star_found == False:
rand_cords = generate_random_coords()
r = rho(rand_coords)
rand_value = generate_random_value()
if rand_value < rand_coords:
return star

# Generate a given amount of stars
def gen_n_stars(n):
for i in range(0, n):
s1.append(gen_star())

number_of_stars_to_gen = ...

# Generate the stars
gen_n_stars(nuber_of_stars_to_gen)



## Forces

So lets start small: We've got two stars whose gravitational pull is acting on each other. The force acting can be calculated using Newton's law of universal gravitation:

$$\vec{F_{12}} = - \frac{Gm_1 m_2}{|r_{12}|^2} \hat r_{12}$$

Symbol Name Value
G Gravitational_constant $\approx 6.674 \cdot 10^{-11} N \cdot kg^{-2} \cdot m^{2}$
$m_1$ Mass of the first object variable
$m_2$ Mass of the second object variable
$r^2$ Distance: $\mid r1 - r2 \mid$ variable
$\hat{r}$ unit vector $\frac{r_2 - r_1}{\mid r_2 - r_1 \mid}$ variable

Now how does this enable us to simulate how galaxies evolve? Well, lets start small, below you can see an image of the forces acting when two stars attract each other.

Now if we want to find out where a star is after one timestep $t$, we need to see how long and how strong a force is acting on the star. So the first thing we need to do, is to calculate the force acting on the star (the overall pull create by alot of other stars). This has to be done for every timestep creating a problem: this is just to much computation that has to be done.

Lets think about the problem for a bit and why it's so important to find a solution to this: Calculating the forces that are acting inbetween 10 stars means calculating how the 100 forces acting inbetween the stars will change their positions. This means that for $n$ stars, the magnitude of $n^2$ stars hast to be calculated and their effects. $n^2$ is alot. Our goal is to get this down to an optimal $n \log(n)$.

This can be done by reducing the amount of stars that are taken into affect by just using the stars in a predefined radius around the star. (I tried this, it was a weird time, lets just skip to the better method). Subdividing the galaxy in cells and calculating the forces acting on a star unsing all the forces occuting from the stars in the cell and the overall forces of the surrounding cells.

So lets get back to the forces: to calculate the new position of a star, we need the value of the force that is acting on in ( $\vec{v}$ ) and the time this force is acting on it ( $t$ ). This is actually simple vector multiplication:

$$F = t \cdot \vec{v}$$

So lets do this with the three stars displayed in the table:

As you can see, we've got three stars with no forces acting on them. They've all got an equal mass. If we use our Equation for calculating the forces inbetween the stars, we get the following results:

## Stability

As you can see, the stars attract each other as expected! At the moment, that's a good sign: we are able to calculate how the stars influence each other, but are still problems to solve... ... one of them is the progrem of stability.
By looking at the stars above, we can see that the stars are all pulling each other towards the area around (2, 2.5). This becomes problematic because we want so simulate a stable system and not a system in which all of the stars get pulled into one point.

The solution is to rotate the galaxy to counteract the force that pulls all stars into the mass center using centrifugal force. There is one problem now: finding out how fast the galaxy can be rotated before it explodes (and how fast it needs to rotate before it implodes).

But before we start creating a stable galaxy, we should find out how to accelerate the calculation and how to display the results:

## Some more code

So what we're doing at the moment is calculating the forces inbetween all the stars in the galaxy. Let's setup some pseudocode:

starsArray = []Stars
starsArray.ImportData("data.csv")
starsArray.CalcAllForces()
starsArray.Draw("out.png")


We're creating an array for storing the stars, adding some stars from a .csv, calculating the forces acting and finally drawing the output.

What does this look like? Something like this:

Looks pretty nice doesn't it? So lets see what we've got: the stars are displayed as a white dot and the force pulling on them is displayed using a line. The line can be seen as a vector showing in what direction the star is pulled. After trying out displaying the force of the vector using its length, I came to the conclusion that it made everything look weird, because of a few stars interacting pretty strongly creating weird lines all over the screen. Thats why the force acting on the stars is diplayed using colors: Red stands for a very big force and blue for a small force. As you can see, the vectors of the stars around the bigger white dots are redder that the ones further away. This means that thre force acting on them is stronger.

## Parallelizing

The image that you can see above containes about 50000 Stars (the image is just an excerpt). Generating it took about 15min utilizing one of the eight threads the Intel® Core™ i7-8550U has to offer. Utilizing go threads, this could be improved nicely:

func main() {
// Define an array storing the stars
starsArray = []Star

// Calculate all the forces in the given galaxy using 8 threads
starsArray.CalcAllForces(8)

...
}

// CalcAllForces calculates all the forces acting inside the Galaxy
// Define a channel for communicating with the workers
channel = make(chan Star, 1000)

// Define how many stars each thread gets to calculate

// Create the given amount of processes
for i = 0; i < threads; i++ {

// Define in what range in the starsArray the thread should calculate the forces
start = i * localRange
end = (i * localRange) + localRange

go calcForces(starsArray, start, end, channel)
}
}


As you can see in the code above, we create an array storing star structs and calling the method that calculates the forces for all the stars in that particular array.

In the CalcAllForces method, we first create a channel that we can use to get the stars from the workers. The channel is buffered, incase of us not being able to process the amount of stars incoming.

We then continue by calculating how many Forces each worker (thread) is going to calculate.

Then, we iterate over the amount of threads staring go-methods that all calculate the forces for a differet slice of stars in the starsArray.

The next task is getting the stars from the workers. This can be done by reading from the channel in the following way:

...

newArray = []Star

for j = 0; j < len(starsArray); j++ {
newStar Star = <- channel
newArray.append(newStar)
}

return newSlice

...


A new array is initialized for storing the newStars that contain information about the force acting on them.

Then, we iterate over the amount of stars (because we should get the amount of stars back from the workers). Because of channels blocking untila you read from them, we can get the new star from the channel, write it to the new array and directly wait for a new star or if there is allready one, write it to the new array.

Because of iterating over the amount of stars, anything that is written below the array only gets called when all the calculations are finished. Because of that, the loop can be seen as a go-thread syncronizer.

The image that is now somewhere up there (the colorful one containing 50000 Star) originally took about 15 min. to calculate using 1 thread. After updating the code, the image can be calculated and generated in around 2 minutes (around 7.5 times faster!!!).

Because of this, we can generate image sequences aka movies / videos in which we can see the motion of the stars.