# I'm trying to simulate moon orbit around earth, my moon is broken.

Then you really need to change the method you are using to integrate the equations of motion. You're basically using Euler's method, which is known to be unstable.
Euler's method is not always unstable; it is so only in some circumstances. The thing to do would be to actually test it and see if it is unstable.

Hurray! The moon works!

Code:
    Public moon_px As Double 'Position.x (m)
Public moon_py As Double 'Position.y
Public moon_pm As Double 'Position Magnitude

Public moon_vx As Double 'Velocity.x (m/s)
Public moon_vy As Double 'Velocity.y
Public moon_vm As Double 'Velocity Magnitude

Public moon_ax As Double
Public moon_ay As Double
Public moon_am As Double

Public angle As Double 'Angular position relative to earth, 0 = "north", 180 = "south" (degrees)

Public gravity As Double
Public earth As Double

Public Timestep As Double '(s)
Public Time As Double

Public Iterations As Double
Public i As Double

Public Degree As Double
Public PI As Double

'---------

Sub DoPrint(thing As Double)

Cells(i + 1, 4).Value = Round(Time / 60 / 60 / 24, 1)
Cells(i + 1, 5).Value = Round(angle * Degree, 1)
Cells(i + 1, 6).Value = Round(moon_pm, 0)
Cells(i + 1, 7).Value = Round(moon_vm, 1)
Cells(i + 1, 8).Value = Round(moon_ax, 3)
Cells(i + 1, 9).Value = Round(moon_ay, 3)
Cells(i + 1, 10).Value = Round(moon_am, 3)
Cells(i + 1, 11).Value = Round(moon_vx, 1)
Cells(i + 1, 12).Value = Round(moon_vy, 1)
Cells(i + 1, 13).Value = Round(moon_px, 0)
Cells(i + 1, 14).Value = Round(moon_py, 0)

Cells(i + 1, 16).Value = angle
Cells(i + 1, 17).Value = Sin(angle)
Cells(i + 1, 18).Value = Cos(angle)

End Sub

'---------

Sub TransLunarInjection()

Degree = 57.2957795
PI = 3.14159265359

moon_px = Cells(2, 2).Value
moon_py = Cells(3, 2).Value
moon_vx = Cells(4, 2).Value
moon_vy = Cells(5, 2).Value
moon_ax = 0
moon_ay = 0
Time = Cells(6, 2).Value
Timestep = Cells(7, 2).Value
Iterations = Cells(8, 2).Value
i = 1

gravity = 6.67384 * 10 ^ -11
earth = 5.97219 * 10 ^ 24

moon_vm = Sqr(moon_vx ^ 2 + moon_vy ^ 2)
moon_pm = Sqr(moon_px ^ 2 + moon_py ^ 2)
moon_am = Sqr(moon_ax ^ 2 + moon_ay ^ 2)

Do While i <= Iterations

If i Mod 1 = 0 Then Call DoPrint(0)

Call Moon_Angle(moon_px, moon_py)

moon_vm = Sqr(moon_vx ^ 2 + moon_vy ^ 2)
moon_pm = Sqr(moon_px ^ 2 + moon_py ^ 2)
moon_am = Sqr(moon_ax ^ 2 + moon_ay ^ 2)

moon_ax = Sin(angle) * ((gravity * earth) / (moon_pm ^ 2)) * -1 * Timestep
moon_ay = Cos(angle) * ((gravity * earth) / (moon_pm ^ 2)) * -1 * Timestep

moon_vx = moon_vx + moon_ax
moon_vy = moon_vy + moon_ay

moon_px = moon_px + moon_vx * Timestep
moon_py = moon_py + moon_vy * Timestep

Time = Time + Timestep

i = i + 0.001

Loop

End Sub

'---------

Sub Moon_Angle(x As Double, y As Double)

If y > 0 And x = 0 Then angle = 0
If y < 0 And x = 0 Then angle = 180 * Radian
If x > 0 And y = 0 Then angle = 90 * Radian
If x < 0 And y = 0 Then angle = 270 * Radian
If x = 0 And y = 0 Then angle = 0
If y > 0 And x > 0 Then angle = Atn(x / y)
If y < 0 And x > 0 Then angle = Atn(x / y) + 180 * Radian
If y < 0 And x < 0 Then angle = Atn(x / y) + 180 * Radian
If y > 0 And x < 0 Then angle = Atn(x / y) + 360 * Radian

End Sub

There were a couple of changes I made thanks to you guys:

- I fixed the radian and degrees issue.
- I fixed my "atan2" Moon_Angle function to properly give angle depending on the moon's position.
- I added a second Timestep multiplyer at my position = position + velocity.
- I Changed my Timestep to 60 seconds.

And it works grand!

I also made sure that my code only prints the results every 1000th loop, as otherwise I ended up with about 40000 lines of data just for one orbit. Now it prints one orbit in 40, give or take a few.

And I think that for my current back-of-the-envelope method of doing this simulation Euler integration will be just fine. The moon makes a perfectly neat orbit and ends up right about where it started. I'm sure that if I made this simulation run for a couple million years something would mess up, but most trans-lunar injections (luckily) don't take millions of years.

I've supplied the datasheet of me simulation for you guys to take a look at, and I thank you all. Your help was greatly appreciated!

#### Attachments

• Moon Orbit Simulation Results Big.pdf
113.9 KB · Views: 122
BvU
Homework Helper

Sub Moon_Angle(x As Double, y As Double)
angle = WorksheetFunction.Atan2(x, y)
end Sub​
and
moon_ax = -Cos(angle) * ((gravity * earth) / (moon_pm ^ 2))
moon_ay = -Sin(angle) * ((gravity * earth) / (moon_pm ^ 2))

moon_vx = moon_vx + moon_ax * Timestep
moon_vy = moon_vy + moon_ay * Timestep

moon_px = moon_px + moon_vx * Timestep
moon_py = moon_py + moon_vy * Timestep​

Works like a dream: zooms around in 28 days in a near-perfect circle with a timestep of 3600 seconds ! Even after 10000 hours, speed is only 0.5% low.

Enjoyed playing with this, even though I just noticed I was a little too late !

BvU
Homework Helper
Hurray! The moon works!

Code:
    Public moon_px As Double 'Position.x (m)

moon_vx = moon_vx + moon_ax
moon_vy = moon_vy + moon_ay

moon_px = moon_px + moon_vx * Timestep
moon_py = moon_py + moon_vy * Timestep
You want to update the vx and vy using the * timestep too !

Gold Member
Is excel this powerful? I though it's something accountants use for their work. -_-

Enjoyed playing with this, even though I just noticed I was a little too late !

Ah well, no worries! Thanks anyway :D

I toyed around with my moon and set the initial position and speed to the moon's perigee and max velocity: http://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html

Ran it like that and I got a very neat eclipse orbit with it's apogee nearly exactly at where our moon's actual apogee is.

Now to start coding my moon ship so I can get a Hohmann transfer underway!

Is excel this powerful? I though it's something accountants use for their work. -_-

It's got a whole programming language under the hood that can be utilized. Of course, it's not the best qua performance to run an astrophysics simulation, but Excel's all we have at school, and it works perfectly fine!

Borek
Mentor
You want to update the vx and vy using the * timestep too !

He did, read the code again.

BvU
Homework Helper
He did, read the code again.
You're right. I had moved the *timestep down to the vx, vy expression to keep dimension of a in lline with naming and also to get the right values for ax, ay, am in the table. In fact, now I start wondering if the expression for px, py improves if 0.5 * timestep is used instead of *timestep...

BvU