Gal_Mond.bas: Galaxy Simulation in 2 D with MOND

Introduction and Purpose

The purpose of this program is to simulate a Galaxy using MOND. The same program also allows you to use Newton's Law.

This program is base around the question: MOND - What is involved
4 Different rotation curves are possible.

The simulation is done using MOND. For an overview of MOND see: "A modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis" by Milgrom, M: 1983ApJ...270..365M

Operation

Input to the program are 12 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
  9. Galaxy Rotation curve ? There are 4 possible rotation curves. 1 means flat. See below for other curves.
  10. Mond ?. Yes means simulation means MOND modification. No means Newton's Law.
  11. Mond Initialization ? This question is only raised if previous question is Yes. No. means that first initialisation cycle does not use MOND.
  12. Universal Constant a0 times 10^8 in cm/sec. Standard value of a0 is 8 * 10-8 cm/sec.
    In the program the value 8 * 10-13 km/sec is used.

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 

200 |xxxxxxxxx                   *               x            **  
    |                     *                    x            *      *
    |                 *                      x             *           *
    |              *                       x              *
    |             *                      x                *
    0             *                    x                  *
      curve 1           curve 2             curve 3           curve 4

                    Galaxy Rotation Curves

Program: GAL_MOND.BAS source

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

This zip file also contains the Visual Basic program: VB Gal MOND.exe For the operation information of this program select: VB Gal Mond operation

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 12 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 : 1
    7. rotation speed v in km/sec: 200
    8. integration methode: 1 (Simpel)
    9. Galaxy Rotation curve: 2
    10. Mond: Yes
    11. Mond Initialisation: Yes
    12. Universal Constant a0 times 10^8 in cm/sec: 8
  2. The program GAL_MOND.BAS is a copy of program GAL_2D.BAS.
    For a general discussion about the other sections See Program "Gal 2D.bas": Galaxy Simulation in 2D
  3. Newton's Law is based around the formula a=G*M/ r^2
    MOND is based around the formula a^2/a0 = G * M/ r^2 or a = SQR(G * M * a0/r^2)
In order to investigate MOND let us first study Newton's Law.
The following figure shows a simple case 
                  <----a----><----a---->
     m2/v2      m1/v1       0        m1/v1     m2/v2
star   4          2         0          1         3
       <----------b---------><---------b--------->

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

           Example 1 Newton's Law, with 2 rings Each with 2 stars.
The above shows 2 linear equations in m1 and m2, which directly can be solved
Now using MOND. Only the ax component has to be used. The factor a0 is temporary removed 
                  <----a----><----a---->
     m2/v2      m1/v1       0        m1/v1     m2/v2
star   4          2         0          1         3
       <----------b---------><---------b--------->

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

         Example 2 MOND, with 2 rings Each with 2 stars.
With MOND we get 2 non lineair equations in m1 and m2. To find a solution is much more complex.

To find a solution goes in three steps:

  1. First we assume that the equations are lineair and are the same as in the case of Newton's Law. With one exeption: the right hand is not equal to acel = v1^2/r but acel^2/a0. This gives an initial value for m1 and m2 (etc). Each of the values v1, v2 etc represent the target rotation curve.
  2. Second we calculate the left side of each equation (ring) and the right side. The left side is called fx and the right side temp. fx/temp is called a correction factor. There is one for each ring.
    When all the correction factions are calculated each mass value is divided by its corresponding correction factor to yield new mass values.
    With fx you can also calculate a speed vx using fx = vx^2/r or vx = sqr(fx*r). vx reflects the actual rotation curve.
  3. Third based on the target rotation curve and the actual rotation curve we calculate an error factor.
    When this error is decreasing step 2 and 3 are repeated. When the error is increasing the following message is displayed:
    Enter E to End C to Exit Loop
    Together with a loop count and error value. When Enter is selected the above calculation is repeated i.e. Step 2 and 3. When C (Continue) is entered the program continues the same as in GAL_2D.BAS
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.

Flat Galxy Rotation curve

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 + G * m0/a^2 = v1^2/a          acel = SQR(G*m1*a0/(2*a)^2) + SQR(G*m0*a0/a^2) = v1^2/a 
 
Example 3 Newton's Law, with 2 stars                  Example 4 MOND, with 2 stars
        and one star in center                               and one star in center
The general solution of example 3 using 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)

The solution of example 4 with 2 stars is:

f(1,1) = G*a0/(2*a)^2
SQR(m1 * f(1,1)) + SQR(G * m0 * a0/a^2) = v1^2/a
SQR(m1 * f(1,1)) = v1^2/a - SQR(G * m0/a^2)
m1 = (v1^2/a - SQR(G * m0/a^2))^2 / f(1,1)
Suppose in Example 4 that m1 = 0
SQR(G * m0 * a0/a^2) = v1^2/a
SQR(G * m0 * a0) = v1^2
v1 = SQR(SQR(G * m0 * a0)
That means in case m1 is zero we get a flat rotation curve
That is also the case when all the mass values in outer rings are zero.

Program Evaluation

The following are tests are all done without any real simulation.
Ater the message:
Perform Simulation, Yes No (E = End)
Select NO.
  1. Execute program with all standard values. Galaxy rotation curve = 2.
    After Mass Calculation phase loop counter = 24. Error = 8,9425
    Observe rotation curve. Observe total mass = 58447645.
  2. Execute program with all standard values. Galaxy rotation curve = 3.
    After Mass Calculation phase loop counter = 45. Error = 0,048
    Observe rotation curve. Observe total mass = 46229844.
    The important lesson is: This is the most accurate curve, which can be simulated using MOND.
  3. Execute program with all standard values. Galaxy rotation curve = 4.
    After Mass Calculation phase loop counter = 22. Error = 58,329
    Observe rotation curve. Observe total mass = 73981419.
    Observe that the highest values of the galaxy rotation curve are 197.1 198.3 198.7 199.0 and 199.1
    The important lesson is: Galaxy rotation type 4 gives a problem See next test.
  4. Execute program the same as before. Galaxy rotation curve = 4 but MOND = Off, That means Newton's Law is used. Observe total mass = 86521296231
    The 6 highest values of the galaxy rotation curve are 193.6 190.0 184.0 176.3 167.7 and 158.7 . That means the speed decreases.
      What is the most important lesson?
    • With MOND it is not possible to simulate a rotation curve which speed decreases.
  5. Execute program with all standard values. Galaxy rotation curve = 1.
    Observe rotation curve. Observe total mass = 1248854978
    The galaxy rotation curve is flat.
    Most important lesson: There is no mass in the disc
  6. Execute program with all standard values. Galaxy rotation curve = 1 (flat) and Rotation speed = 140
    Observe rotation curve has the values 140.0 200.0 200.0 200.0 etc.
    Observe total mass = 299850080
  7. Execute program with all standard values (with MOND) except Mass of Black hole = 15000000 and Galaxy rotation curve = 1 (flat).
    After Mass Calculation phase loop counter = 1000. Error = 1,227
    Observe rotation curve has the values 200.2 200.3 200.4 200.4 etc. Observe total mass = 15000003873
    The mportant lesson is: By introducing a large BH you can get a flat galaxy rotation curve
  8. Execute program with all standard values except Mass except Galaxy rotation curve = 1 (flat) and MOND select N. That means use Newton's
    Observe total mass = 130885015144
      What can you learn from the tests 5, 6, 7 and 8 ?
    • First that you can simulate flat galaxy rotation curves with MOND but that the mass of such a galaxy with MOND is a factor 100 smaller than with Newton's Law.
    • Second the best solution with MOND is a black hole at the center but still a the mass a roughly a factor 10 less.
  9. Execute program with all standard values except Mass except Galaxy rotation curve = 1 (flat). Select MOND but when you get the message "MOND intialisation" select NO
    you will get message "Loop 1 Error 2916.604" and the message "Enter E to End C to Exit Loop". hit C (Continue) immediate.
    Observe total mass = 130885015144. Observe rotation curve.
    Observe starting value of 315.7, intermediate value of 1013.5 and final value of 1576.8
      What can you learn from test 8 and 9 ?
    • First, the main point is that the total mass and the mass distribution is identical but the rotation curve not With Newton's Law the curve is flat but with MOND the speed increases almost lineair starting from 300 to roughly 1600
    • Secondly that the rotation curves, assuming the same baryonic mass distribution of a simulation based on Newton's Law compared to MOND are totally different.
    • Third As already mentioned with MOND at larger distances the speed increases or stays constant, but does not decrease.
  10. Execute program with all standard values. Galaxy rotation curve = 2. Except Universal constant a0 = 0.8
    When you get the message: "Enter E to End C to Exit Loop" hit Enter.
    Do this 10 times until you get message "Loop 61 Error 2.841". Observe that error is decreasing.
    After 10 times hit C to Continue. Observe rotation curve. Observe total mass = 599452127.
    Repeat same test with a0 = 0.08. Observe total mass = 5994521270.
      What can you learn from the test 9 and 10 ?
    • In this case we keep the rotation curve the same.
      That means the above formula applies: v1 = SQR(SQR(G * m0 * a0)
      With v1 being constant when you decrease a0 with a factor 10, m0 increases with a factor 10. And vice versa.
    • That by changing the Universal Constant a0 you can give the galaxy any mass value you like.
  11. Now let us study what happens when you keep the mass constant.
    Execute the same test as #9. Mond Initialisation = NO. Universal constant = 0.08
    you will get message "Loop 1 Error 564.370". Hit C (Continue) immediate.
    Observe total mass = 130885015144. Observe rotation curve.
    Observe starting value of 99.8, intermediate value of 320.5 and final value of 498.6
  12. Do the same test again. Universal Constant = 0.0008 you will get message "Loop 1 Error 331.256". Hit C (Continue) immediate.
    Observe total mass = 130885015144. Observe rotation curve.
    Observe starting value of 31,6, intermediate value of 101.4 and final value of 157.7
      What can you learn from the tests 9, 11, 12?
    • In those three tests we keep the total baryonic mass constant. MOND is tetsted by modifying a0
    • With Newtons'Law we get a flat curve with speed 200.0
    • With MOND we get a curve which lies below this curve
    • IMO what this means is to use MOND in order to explain dark matter is highly speculative.
The following are two tests with a real simulation.
Ater the message:
Perform Simulation, Yes No (E = End)
Select Yes.
  1. Execute program with all standard values except # of rings is 5 and number of stars = 20. MOND = Yes. When you get the message: "Enter E to End C to Exit Loop" hit Enter.
    Do this 10 times until you get message "Loop 86 Error 11.027". Observe that error is decreasing.
    After 10 times hit C to Continue. Observe total mass = 54078319
    Select Yes
    The Sun is at ring 2. Number of stars = 100. The first signs that chaos starts is at angle 600. Chaos becomes visible at 700 but stars still move in circles around 1000. No stars are ejected. Hit "N" to terminate simulation.
  2. Execute program with all standard values except # of rings is 5 and number of stars = 20. MOND = No. Observe total mass = 154507709241. This is a factor 300 larger as with MOND! Same Rotation Curve!
    Select Yes
    The Sun is at ring 2. Number of stars = 100. The first signs that chaos starts is at angle 200. Chaos becomes visible at 210. Real chaos at 230. Stars are ejected. Hit "N" to terminate simulation.
  3. Repeat both tests but now also observe the variable a. A is the maximum value of acceleration.
    In test 1 the initial value is 0.015. At angle 600 this becomes 0.1. Around 700 0.3 and final values are around 1.
    In test 2 the initial value is 0.47. At angle 200 this becomes 0.55 and final values are around 100.
      What is the main lesson ?
    • That simulations solely with MOND are more stable than with Newton's Law.
    • That a more realistic simulation should be based on both Newton's Law and MOND. But how do you define the regio where both apply?
    • That when chaos starts with MOND no stars are ejected.
Created: 13 March 2007
Modified 14 August 2015

Back to my home page Contents of This Document