Over the previous year, I learnt about physics simulation and a lot about physics in general, mainly from brilliant. Thanks to those who helped me a lot along the way, Neeraj Anand Badgujar (Talulah Riley), Steven Chase, and Karan Chatrath.
First, we need data for each position of the bodies; you can do this by assigning random initial positions to the bodies within a certain range.
Build arrays that contain mass data of these bodies, and components of the position values, and and components of the velocities.
Now, here's how this brute-force CPU simulation works:
First, let's just try to figure out the dynamics of one body, for a small timestep.
Forces can be computed by multiplying the unit vectors to each other particle and the magnitudes of the gravitational force.
Essentially, 3 things should be computed each time for the bodies to find the total force exerted onto only 1 of these bodies.
The distance between the single body and each other body, the unit vectors, and finally gravitational force:
To do this as a program, loop through each of these bodies, and keep a total of the force which you keep updating until you reach the final th body.
Then integrate using Euler integration.
This should be done to all the masses. Hence it will result in a nested loop. However, looping from 1 to the th body results in a problem. When you start with the first body, you could be calculating the distance from the first body to the first body, which is zero. Then if you divide by the distance for the unit vector, this would result in division by zero.
So, you should only iterate the algorithm when the indexes of the position data and are not equal, ie. so you skip the the th body when doing force computation.
Another consideration: what will happen if both bodies collide perfectly; ie. their centers touch? This would result in division by zero as well in the force computation; we cannot allow that to happen.
In my code I just let each body pass through one and other when they are a certain distance away from each other; this is done by setting the force to zero when they are a small distance away from each other, therefore they just move as per the previous values of their velocity.
However, once they are no longer in that small range, gravitational force again takes precedence.
If you want you can code in some resistive forces, like electrostatic force where all the charges are positive or all are negative, you can just calculate it along with the gravitational force using Coulomb's law. You can also add collision mechanics to it to prevent this error.
Also, when you want to animate each body, it will be a little computationally expensive. I found that using the Matlab rectangle() function to draw a filled circle less expensive than viscircles(). All you have to do is draw the rectangle with curvature and add a random colour to it. I personally added random colour using a colour matrix for random RGB values with row vectors for each body.
Like most of my simulations, I put the dynamics data into arrays which correspond with the time frame of the system and animate separately to the simulation.
However, since there are bodies, we should have 2 matrices where is the number of frames the animation has, one for each component in and .
The computation does run on the CPU, so wait a few seconds before the animation happens.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122