N-Body - Moon won't stay in orbit

[neoaikon]

A long time ago, I wrote a program that takes bunches of particles and applies Newtonian gravity to each so they interact. Lately I've rewritten this program, in an effort to simulate planetary motion. However I've ran into a bit of a snag, I can get a Sun-Earth system to behave properly, and I can get a Earth-Moon system to behave properly, but when a Sun-Earth-Moon system is created, the Earth will orbit fine, but the moon will simply follow the Earth around, never orbiting around it, and the distance between the Earth-Moon will slowly increase. I have attached 3 images (Absolutely not to scale):

The first shows the simulation right after start, the Gray dot is the moon, slightly under it but somewhat visible is the Earth in Blue, and the Yellow orb is...well I'll let you guess ;)

The second shows the moon clearly trailing behind the earth. And the third for good measure is after the moon has reached about 1AU from the sun.

To calculate the initial velocities I use the formula v = gc * (m / d), I apply this velocity on a vector perpendicular to the main orbiting body (Sun for Earth, Earth for Moon). For the moon I'm calculating the velocity to orbit the earth, and then adding the velocity to orbit the sun to make it relative to the sun.

I find the Gravitational Acceleration using the formula a=gc * (m / d^2), I then find a normalized 2d vector to the gravitating body. I multiply this acceleration by the components of the 2d vector. and then add it to the XY position. For the moon here is the algorithm:

- Calculate Vector to Earth
- Calculate Vector to Sun
- Calculate the Acceleration from Earth
- Calculate the Acceleration from Sun
- Moon's X = X + (SunAcc * -(SunDirX)) + (EarthAcc * -(EarthDirX))
- Moon's Y = Y + (SunAcc * -(SunDirY)) + (EarthAcc * -(EarthDirY))

Whew, that's a lot to type. Anyways, as you can see in the attached images, the moon doesn't exhibit the expected motion. It seems that the way I am adding the acceleration to the moon object is in error, but I have checked and recheck my formulas, masses, distances, vectors and acceleration for errors and cannot seem to spot it. I appreciate any and all help! This has me going in circles it seems.

 

[neoaikon]

Looking back on some of my past notes I realized that the forces need to be added together before being multiplied into the vectors. The final moon routine looks like this:

- Calculate Vector to Earth
- Calculate Vector to Sun
- Calculate the Acceleration from Earth
- Calculate the Acceleration from Sun
- Moon's X = X + (TotalAcc * -(SunDirX)) + (TotalAcc * -(EarthDirX))
- Moon's Y = Y + (TotalAcc * -(SunDirY)) + (TotalAcc * -(EarthDirY))

Unless this is incorrect, this improves the behavior of the moon. It now will stay tethered to AND orbit around it, but slowly the moon still drifts away . Because of the vast improvement in behavior, I think it may be how I'm starting the simulation.

I start with the sun in the middle, the Earth is directly to the left of it with its vector pointed downward and perpendicular to the sun. The moon is to the left of this, with its vector pointed downward and perpendicular to the earth. I then calculate the initial velocities, multiplied by the vector, to set the planets up. So the Earth goes counter-clockwise around the sun, and the moon counter-clockwise around the earth.

The moon is the only object receiving gravity from 2 objects, the Earth is only effected by the sun, and the sun doesn't move at all. And like I said, I checked my starting velocities, masses, gravitational constant, vectors, the works, so it has to be in how I'm handling the acceleration or how I'm positioning everything initially.

 

[atyy]

Do you get the same result if you make the time step say 10X or 100X smaller? If you don't, then the inaccuracy is in the numerical accuracy of the simulation. If you do, then maybe the equations you used are incorrect. (There's also the possibility that the system is genuinely unstable, but I'd think about that later. BTW, your attachments haven't been approved, so I haven't seen them.)

 

[neoaikon]

I do get the same results actually. The system should be stable, according to all other simulations I've seen, since I'm using real world values for these things. I wish I had the original source to peek at, all I have is my notes.

I'm pretty sure the equations are correct, I've used them before and they're in my handy notebook of things like that. Whatever is happening is within how I am applying the forces to the moons direction vector, I think. The difference between the moon and the Earth is there are 2 direction vectors, if I remove the sun, and keep the moon relative to the earth, it will orbit fine, however that's not what I want since I love the beauty of dynamic systems. So I'm sure the force is just being applied in the wrong direction, Since I've noticed the orbit of the moon becomes very eccentric (now that its orbiting at all)

I'll be happy when I figure this all out :uhh:

 

[D H]

Looking back on some of my past notes I realized that the forces need to be added together before being multiplied into the vectors. The final moon routine looks like this:

- Calculate Vector to Earth
- Calculate Vector to Sun
- Calculate the Acceleration from Earth
- Calculate the Acceleration from Sun
- Moon's X = X + (TotalAcc * -(SunDirX)) + (TotalAcc * -(EarthDirX))
- Moon's Y = Y + (TotalAcc * -(SunDirY)) + (TotalAcc * -(EarthDirY))

This isn't valid. Look at it in terms of units. X and Y are positions: they have units of length. Acceleration is length/time squared. You can't add position and acceleration any more than you can add any other pair of dissimilar quantities such as distance and mass. What is 1 meter plus 1 kilogram?

What you need to do is use the velocities. Acceleration changes the velocity as time changes, and velocity changes the position as time changes. Some nomenclature first:

$$\aligned \boldsymbol{x}_m(t)\quad&\text{is the time-varying position vector of the moon} \\ \boldsymbol{v}_m(t)\quad&\text{is the time-varying velocity vector of the moon} \\ \boldsymbol{a}_m(t)\quad&\text{is the time-varying vector acceleration} \\ r_{m,s}\quad&\text{is the time-varying distance between the moon and the sun} \\ \hat{\boldsymbol{r}}_{m,s}\quad&\text{is the time-varying unit vector from the moon to the sun} \\ r_{m,e}\quad&\text{is the time-varying distance between the moon and the earth} \\ \hat{\boldsymbol{r}}_{m,e}\quad&\text{is the time-varying unit vector from the moon to the earth} \aligned$$

The Earth and Sun have similar states.

Here is a very simple state propagator: First calculate the accelerations of the Moon, the Earth, and the Sun toward each other. The Moon's acceleration toward the Earth and Moon is

$$\boldsymbol{a}_m(t) &= \frac{GM_s}{r_{m,s}^2(t)}\hat{\boldsymbol{r}}_{m,s} + \frac{GM_e}{r_{m,e}^2(t)}\hat{\boldsymbol{r}}_{m,e}$$

The accelerations of the Earth and Sun are similar. Now update the states to your next simulation time step:

$$\aligned \boldsymbol{x}_m(t+\Delta t) &= \boldsymbol{x}_m(t) + \Delta t\,\boldsymbol{v}_m(t) \\ \boldsymbol{v}_m(t+\Delta t) &= \boldsymbol{v}_m(t) + \Delta t\,\boldsymbol{a}_m(t) \endaligned$$

The approach outlined above is too simple. It is called Euler integration and is not numerical stable. Just as simple, but a lot more stable, is to calculate the new velocities first:

$$\aligned \boldsymbol{v}_m(t+\Delta t) &= \boldsymbol{v}_m(t) + \Delta t\,\boldsymbol{a}_m(t) \\ \boldsymbol{x}_m(t+\Delta t) &= \boldsymbol{x}_m(t) + \Delta t\,\boldsymbol{v}_m(t+\Delta t) \endaligned$$

This simple change will make your simulation a lot more accurate. It is called Euler-Cromer integration by some, symplectic Euler integration by others. Welcome to the world of numerical integration! Euler and symplectic Euler are just the tip of the iceberg. There are much better schemes that are also more complex, some considerably more complex.

 

[neoaikon]

Wow! Thanks for the info D H. I was actually hoping to draw your attention, because a previous post of yours helped me out, and that's what made me choose here to post these types of questions. :). I'll add the changes you specified and see if that makes it more stable.

 

[neoaikon]

Well I actually am using velocity, but the terminology within my shortcode is wrong.

- Calculate Direction Vector to Earth
- Calculate Direction Vector to Sun
- Calculate the Acceleration Vector from Earth and Sun
- Calculate the Velocity Vector from the Earth and Sun
- Moon's X = X + (SunVel.X * -(SunDirX)) + (SunVel.X * -(EarthDirX))
- Moon's Y = Y + (SunVel.Y * -(SunDirY)) + (SunVel.Y * -(EarthDirY))

That is more accurate. If anyone wants to look at my code I can upload it to my webserver, its in C#.

 

[edguy99]

Well I actually am using velocity, but the terminology within my shortcode is wrong.

- Calculate Direction Vector to Earth
- Calculate Direction Vector to Sun
- Calculate the Acceleration Vector from Earth and Sun
- Calculate the Velocity Vector from the Earth and Sun
- Moon's X = X + (SunVel.X * -(SunDirX)) + (SunVel.X * -(EarthDirX))
- Moon's Y = Y + (SunVel.Y * -(SunDirY)) + (SunVel.Y * -(EarthDirY))

That is more accurate. If anyone wants to look at my code I can upload it to my webserver, its in C#.

Be careful what you use for starting values of velocities. Try tilting the starting velocity angle in a bit as using a tangent to the circle for starting will slightly overstate your 'outward' velocity and this will slowly compound.

Hope this is useful.

 

[neoaikon]

Shouldn't the results of my formula provide good starting velocities? The orbits I get from that seem to fall within range, Is there anyway I can find out how much I should move the starting angles inward? or the vector needed for said velocity?

 

[edguy99]

Shouldn't the results of my formula provide good starting velocities? The orbits I get from that seem to fall within range, Is there anyway I can find out how much I should move the starting angles inward? or the vector needed for said velocity?

Depends on the time slice you are using for calculations (do you calculate once per minute, once per hour, or once per day or ?). To get the starting velocity vector, visualize where the moon was at time 0 minus 1, ie. it was not directly above where it is at time 0 (as is implied if you use the tangent and point straight down). The smaller the time slice you use, the smaller this error is, but you can get correct orbits at almost any time slice if you fix the starting vector - good luck.

 

[neoaikon]

I have came to the conclusion that Euler Integration will not work for me in this case, as the low compound percision of the measurements seems to affect the planetary system in a bad way.

I am now using Velocity Verlet Integration, as I've read good things about it. However I'm not familiar with the fancy equations, so I'm having a bit of trouble learning how to implement it properly. Currently I have it working with the following algorithm:

-Set Initial Acceleration and Position
-EarthX = (2 * EarthX - OldEarthX) + (0.5 * OldEarthXAcceleration * (dt * dt))
-EarthXAcceleration = (GC * SunMass) / (EarthSunDistance * EarthSunDistance)
-EarthXVelocity = 0.5 * (EarthXAcceleration + OldEarthXAcceleration) * dt
-OldEarthXAcceleration = EarthXAcceleration

I do the same thing for the Y coordinates, but I'm going to be lazy and only type out X. I came to this algorithm reading "Velocity Verlet" section at http://wiki.vdrift.net/Numerical_Integration#Velocity_Verlet"

Is my implementation of this method correct? I'm hoping the stability this method provides should allow me to get my moon working, but while it does give orbital motion in my sim, I don't know if it's correctly done.

Additionally, with this method I can no longer set my starting velocity with v = gc * (m / d) with this method, and I'm unsure of how I should set this value.

Thank you all for your help thus far! Many thanks.

EDIT: I have to divide my initial velocity by 38,677 to get circular motion, and I'm really scratching my head as too why...if anyone wants to look at my code, you can find it at http://aikonland.hyperphp.com/csMain.txt" .

 

[neoaikon]

I have determined my Velocity Verlet implementation to be correct, and the problem with the initial velocities. for an Earth like mass the solution for the most circular orbit was 38,677, and the solution for a Mars like mass was 38,575. Although I brute forced both of those values by trying many different ones till I got the orbits circular. I just need to know how to set up the initial velocity properly.

On a side note I am in awe at how stable velocity verlet is. t=5,000,000 with a DT of .1 and my simulation is still stable and in their original orbits with a Mars-Earth-Sun system. I plan on adding mercury and venus just for kicks once I brute force their values, or find out how to set it up correctly.

 

[neoaikon]

Ha! I figured out the issue, and a few mistakes in my implementation of VV.

- Set Initial Velocity and Position
- EarthX += (EarthVelX * dt) + (0.5 * OldEarthXAcceleration * (dt * dt))
- Get EarthSunDirection Vector
- EarthXAcceleration = ((GC * SunMass) / (EarthSunDistance * EarthSunDistance)) * EarthSunDirectionX

- EarthXVelocity += 0.5 * (EarthXAcceleration + OldEarthXAcceleration) * dt
- OldEarthXAcceleration = EarthXAcceleration

As you can see there is a differece in how I was updating the X position. I included the Direction vector which I was doing but didn't put into the shortcode. and I had = instead of += on the velocity vector too. Now that those mistakes are eliminated, the true reason my initial velocity was getting set incorrectly wasn't in the setting of it, but in the update position formula

- EarthX += (EarthVelX * dt) + (0.5 * OldEarthXAcceleration * (dt * dt))

Simply, I didn't need to add the OldEarthXAcceleration into the EarthX! When I really thought about it, and worked out the math, it only takes 20000 steps for it to make the Earth crash into the sun at that rate. when I removed it, the heaven's aligned, and god smiled.

Now at least I can try to get back to the meat of this question. and see if the moon will freaking orbit properly with VV! On a side note I have a question! Does velocity vector simply average the acceleration across timesteps? looking at the code it that looks like the case.