Program "Gal 2D.bas": Galaxy Simulation in 2 D

Introduction and Purpose

The purpose of this program is to simulate a Galaxy in 2D with Quick Basic.
This program is base around the question: Stable Galaxy - What is involved

The galaxy rotation curve is flat.

Operation

Input to the program are 8 parameters:
  1. # of stars per ring. Minimum number is 2. Maximum Number is 50.
  2. # of rings with stars. Maximum # is 50.
  3. # of steps per revolution. Minimum # is 1000.
  4. Mass of Black hole in the centre in increments of 1000 sun masses. 0 = No Black hole
  5. Extended disc ? (yes or no) ? No. means normal disc size i.e. 60000 light years
  6. Extended disc scale ? This question is only raised when the previous question is answered with Yes. 1 implies no extended disc. Standard value is 10 i.e. 600000 light years
  7. Rotation speed v in km/sec. 100 is Minimum. 400 is Maximum
  8. Integration method 1 = simple 2 = Runga Kutta (K1,K2,K3,K4).
    For details see: http://en.wikipedia.org/wiki/Runge-Kutta
    Or go to http://www.analyzemath.com/ and goto Section 5 "Calculus Tutorials": and select Differential Equations - Runga Kutta Method

The display shows:
  1. The position of all the stars each time after 10 steps.
  2. In the left top corner the number of the ring that harbors the sun.
  3. In the left top corner besides that number the angle that the sun has moved. A value of 360 means that the sun has moved one revolution.
  4. The distance of at least 1 star in each ring
  5. The speed of that same star

Program: GAL_2D.BAS source

In order to download the source select:GAL_2D.BAS
To see the program listing select:GAL_2D.HTM
Execution file select: GAL_2D.zip and: brun45.exe

The same program "Gal_2D" is also written in Visual Basic 5.0
For the operation information of this program select: Visual Basic program "VB Gal 2D" - Description and operation
For the executable of the program select: VB Gal 2D.zip

Technical Data

The program consists of 2 parts: A calculation part and a simulation part.
The calculation part consits of 5 sections.
  1. In the first section the 8 parameters are entered. The standard values are:
    1. # of Rings: 10.
    2. # of stars per ring: 50
    3. # of steps per revolution: 1000
    4. Mass of Black hole in the center: 0
    5. Extended disc halo : No
    6. Scale factor of extended disc : 10
    7. rotation speed v in km/sec: 200
    8. integration methode: 2 (Runga Kutta)
  2. In the second section
    1. distance between the rings is calculated.
      First the distance between ring 2 and ring 1 is 2. The distance between ring 3 and ring 2 is 3. etc.
      Secondly the distance is scaled based on a maximum distance of 60000 light years.
    2. The angle between each star in each ring is identical. The angle(n) of star n in a ring of nstars is equal to n * 2 * pi / nstar
    3. The position (x,y) for each of the stars is calculated. For ring i: x(n) = dist(i) * cos(angle(n)). y(n)= dist(i) * sin(angle(n)).
  3. In the third section the speed of the stars per ring is calculated.
    The (absolute) speed of all the stars in each ring is identical: For ring i : vvr(i) = vvr(1)
    The initial value for the total mass of the stars m(i) in each ring is set equal to 1.
  4. In the fourth section the acceleration for each of nring*nstar stars is calculated and stored in the array f(iv,i)
    This calculation is performed in three loops, with the variables: iv, i and j.
    The sketch below shows this calculation for # of rings = nring = 3 and # of stars/ring nstar = 8
                        19
 

        17          11           18

             12             10 
                     3       
                 4       2    

   21   13    5      0     1     9    17 

                6        8                
                     7        
            14               16  

        22           15            24

                     
                     23 
   
    Example 1 of 3 rings Each with 4 Stars
    one star in each ring (i.e. the stars numbered 1,9 and 17). This is the source star. Only one star has to be calculated because the acceleration of the stars in each ring is identical. The acceleration of star 2 until 8 is identical as star 1 etc
    First the variable iv is calcuted from 1 to nring. This is done to select the source star.
    Secondly the variable i is caluculated from 1 to nring. This is done to select the destination ring. The total acceleration factor in both x and y direction is set to zero.
    Third the variable j is calculated from 1 to nstar. This is done to select the destination star. Using the previous calculated x and y values of the source and destination stars we can now calculate the distance between those 2 stars, the acceleration towards the destination star, (G/r^2) and the acceleration in both the x and the y direction. Those two values are used to update the total acceleration factors.
    When this is done for all the stars it can easily be seen that the total acceleration factor in the y direction is zero. Only the total acceleration factor in the x direction is non zero. That value is stored in the array f(iv,i).
    This is repeated, in order to select a new destination ring, for all the variables i until nring.
    This is repeated, in order to select a new source star, for all the variables iv until nring
    At this point the whole aray f(iv,i) is calculated. 
                      <----a----><----a---->
     m2/v2      m1/v1       0        m1/v1     m2/v2
star   4          2         0          1         3
       <----------b---------><---------b--------->

    G/(2a)^2              * m1   G/(a+b)^2 - G/(b-a)^2 * m2 = v1^2/a
    G/(b-a)^2 + G/(a+b)^2 * m1   G/(2b)^2              * m2 = v2^2/b

          f(1,1)          * m1        f(1,2)           * m2 = v1^2/a
          f(2,1)          * m1        f(2,2)           * m2 = v2^2/b

        
        Example 2 with 2 rings Each with 2 stars.
  5. In the fifth section the mass for each of the stars (per ring) is calculated.
    1. The following sketch shows the situation when one ring is considered. 
                <----a----><----a---->                  <----a----><----a---->
        m1/v1       0        m1/v1              m1/v1       m0       m1/v1  
star      2         0          1                  2          0         1

      acel = G * m1/(2*a)^2 = v1^2/a        acel = G * m1/(2*a)^2 + G * m0/a^2 = v1^2/a 
 
    Example 3 with 1 ring with 2 stars.    Example 4 with 2 stars and one star in center
      The general solution Newton's Law is:
      m1 * f(1,1) + G * m0/a^2 = v1^2/a
      m1 * f(1,1) = v1^2/a - G * m0/a^2
      m1 = (v1^2/a - G * m0/a^2) / f(1,1)
    2. In the case of n rings we have to solve n equations with n variables i.e. the n masses of the individual stars of the n rings. THis is done by using the subroutine SOLVEEQUATIONS. 
                                 <----a----><----a---->
    m3/v3     m2/v2      m1/v1      m0        m1/v1     m2/v2    m3/v3
star  6         4          2         0          1         3        5
                <----------b---------><---------b--------->
      <--------------c---------------><---------------c------------> 

               G/(2a)^2 * m1   G/(a+b)^2 - G/(b-a)^2 * m2   G/(a+c)^2 - G/(c-a)^2 * m3 = v1^2/a - m0*G/a^2
  G/(b-a)^2 + G/(a+b)^2 * m1                G/(2b)^2 * m2   G/(b+c)^2 - G/(c-b)^2 * m3 = v2^2/b - m0*G/b^2 
  G/(c-a)^2 + G/(c+a)^2 * m1   G/(c+b)^2 + G/(c-b)^2 * m2                G/(2c)^2 * m3 = v3^2/c - m0*G/c^2

        f(1,1)          * m1        f(1,2)           * m2            f(1,3)       * m3 = v1^2/a - m0*G/a^2
        f(2,1)          * m1        f(2,2)           * m2            f(2,3)       * m3 = v2^2/b - m0*G/b^2
        f(3,1)          * m1        f(3,2)           * m2            f(3,3)       * m3 = v3^2/c - m0*G/c^2
      

        Example 5 with 3 rings Each with 2 stars and star in centre = m0
      The above sketch shows the situation when 3 rings are involved. In that case 3 masses have to be calculated.
      The total mass of each ring i is nstar times m(i).
At this point of the program, summary data is displayed.
  1. the total acceleration value ax and ay. Observe that ay is small compared to ax.
  2. For each ring: the speed, the total mass and the density.
  3. The ring that harbors the sun.
  4. For the whole disc the total mass.
  5. The revolution time T of the sun and the value of dt. dt = T / # of steps per revolution.
The next question is displayed:
Perform Simulation, Yes No (E = End)
If you enter Y then the simulation is done using Newton's Law.
Solution is calculated in a simple manner (method = 1) or using Runga Kutta method (method = 2)

To stop the simulation: Enter Escape, E to End or N to continue with a New simulation.

Program Evaluation

  1. Execute the program and Select # of Rings is 2, # of stars = 25. Keep all the other parameters standard. Just hit Enter. When the simulation starts you will see 2 rings showing each 25 stars, which nicely move in a circle until the angle is roughly 250 degrees. Then chaos starts.The sun is ring 1
  2. Execute the program and Select # of Rings is 5, # of stars = 20. Keep all the other parameters standard. Just hit Enter. When the simulation start you will get 5 rings showing each 20 stars, which nicely move in a circle until the angle is roughly 150 degrees. Then chaos starts. The sun is ring 2
  3. Try 2 rings with 6 stars. Then sun is ring 1. Chaos starts at roughly 2600 degrees.
  4. Try 2 rings with 4 stars. Chaos starts at roughly 80 degrees.
  5. Try 2 rings with 2 stars. Chaos starts at roughly 10 degrees.
  6. Try 1 rings with 4 stars. Chaos starts at roughly 2400 degrees.
What the simulations show in general that the program works; you get stable galaxies with flat rotation curves.
However at a certain moment you get chaos. There are 2 reasons:
  1. There is a problem with the mathematics involved for systems with 2 rings and a small number of stars
  2. Accuracy.
  • In order to understand the first reason you have to study example 2. Example 2 shows the initial condition of 4 stars all in one line ( 2 rings with each 2 stars) The masses of the stars are calculated based on this configuration. As a result only the x component is taken into account. The y componenten are zero.
    This is also valid if there are more rings and more stars in each ring.
    The problem starts when the stars starts to move. The speed v of all the stars is the same but the angle-speed not. The stars near the centre move over a much larger angle compared with the stars near the rim. The result is that the stars are not more alligned. Example 6 shows that situation. 
                                                   m2/v2
                                                 C

                                       b
                                               b*sin(g) 
                  <----a---->
                m1/v1       0        m1/v1        
      F           B         0          A         E
                             <----a---->

                  b
              
     m2/v2      
      D   
        
        Example 6 with 2 rings Each with 2 stars.

     star 1 = A  star 2 = B  star 3 = C and star 4 = D
     Distance between star A (and B) and centre 0 is a 
     Distance between star C (and D) and centre 0 is b 
     angle C0E = g   angle F0D = g   
     CE = FD = b*SIN(g)   AE = b*COS(g) - a   FA = b*COS(g) + a

     Force of A towards C: FAC = G/ (AE^2 + EC^2) = G/((b*COS(g) - a)^2 + b*SIN(g)^2)
     X component = FAC * AE / AC = FAC * (b*COS(g) - a) / SQR(((b*COS(g) - a)^2 + b*SIN(g)^2)
     Y component = FAC * CE / AC = FAC * b*SIN(g) / SQR(((b*COS(g) - a)^2 + b*SIN(g)^2)

     Force of A towards D: FAD = G/ (AF^2 + FD^2) = G/((b*COS(g) + a)^2 + b*SIN(g)^2)
     X component = FAD * FA / AD = FAD * (b*COS(g) + a) / SQR(((b*COS(g) + a)^2 + b*SIN(g)^2)
     Y component = FAD * FD / AD = FAD * b*SIN(g) / SQR(((b*COS(g) + a)^2 + b*SIN(g)^2)
    For the matrix component f(1,1) (i.e. the force A towards B) there is no difference.
    For the matrix component f(2,1) (i.e. the force A towards C and D) there is a difference.
    First the two y components are non zero and not the same i.e. there is a small force perpendicular to the line "FB0AE". Secondly the values of the x components are different. The result is that the calculated mass values will be different. The good news is that this is only a problem when there are a small number of stars involved.

    To test this situation perform the program and select # of rings = 2 and # of stars is 4. When the message "Perform Simulation, Yes No (E = End)" Enter R for repeat. Do this again and again. Each time observe the mass values for each ring and the values for ax and ay. You will see that the mass values are different and that ay has almost the same value as ax (except for the first time)
    Perform the same test but now for 2 rings and 40 stars. Observe that ay is small compared to ax.

  • The second reason: accuracy, comes into play when the number of rings is 1. In that case all the stars should move exactly with the same angle forward and staying at the same distance from the origin O. This does not happen. There is a small error and After each calculation cycle of all the stars, this error becomes larger and larger until the system becomes unstable.
    • For example: Try 5 rings 50 stars (methode 2, # of steps per revolution 1000). This configuration of 250 stars becomes unstable at angle 71. Because the the system consits of 5 rings the distance of the first 50 stars should be identical and equal to exactly 5800. This is not the case. There is a small error. You have to modify the program to make this visible. What happens there is a small error which slowly increases and the distance is also not the same.
    • Try the same but now with # steps per revolution 2000. Again you get a small error but this error is initially smaller than with the first test. This seems better, but slowly the error increases and becomes the same as with the first test resulting in both cases chaos around angle 71.
    • Try the same (5 rings, 50 stars, # of steps per revolution 1000) but now with with method 1. What you get after each cylce is that the distance changes, but the distances are the same for each ring. Not exactly the same there is a small error between each (with in the same ring) and this error increases after each cycle until again you get chaos angle 71.
modified: 26 January 1999
modified: 23 February 2007
modified: 21 July 2013
modified: 20 August 2015
updated : 19 Februari 2016

Back to my home page Contents of This Document