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

1. May 26, 2014

### xzbobzx

1. The problem statement, all variables and given/known data

Alright, this is probably going to be a long one.

For a school project I am trying to simulate (in 2D) our moon's orbit around the earth, to eventually simulate a spaceship doing a hohmann transfer between LEO and LLO whilst being affected by both the earth's gravity and the moon's gravity.

I am not trying to simulate/solve a three body problem. And the simulation is more of a back of the envelope simulation than a true earth-moon system simulation, but it will serve it's purpose:
The earth is not affected by the moon.
The moon is only affected by the earth's gravitational pull and its initial velocity.

The given variables/data/setup:
- I have a coordinate system X and Y, Y goes up, X goes right. The (static) earth is located at the centre.
- Our moon is located at Y = 385000000 m (the moon's mean orbital distance is 385000 km)
- Moon initial velocity is X = 1023 m/s (the moon's mean orbital velocity is 1.023 km/s)
- Moon mass is irrelevant for now as I'm doing it all based on acceleration, and because the earth is static.
- Earth mass is 5.97219*10^24 kg.
- A Timestep of 1 second.

2. Relevant equations

$g = -\frac{GM}{r^2}$

$g =$ gravitational acceleration
$G =$ gravitational constant
$M =$ earth mass
$r =$ orbit radius

3. The attempt at a solution

I set up a Microsoft excel project and began coding an macro do to all my heavy lifting:

Code (Text):

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

'---------

Sub TransLunarInjection()

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_ax = 0
'moon_ay = 0

Do While i <= Iterations

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_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
moon_py = moon_py + moon_vy

Cells(i + 1, 4).Value = Time
Cells(i + 1, 5).Value = angle
Cells(i + 1, 6).Value = moon_pm
Cells(i + 1, 7).Value = moon_vm
Cells(i + 1, 8).Value = moon_ax
Cells(i + 1, 9).Value = moon_ay
Cells(i + 1, 10).Value = Sqr(moon_ax ^ 2 + moon_ay ^ 2)
Cells(i + 1, 11).Value = moon_vx
Cells(i + 1, 12).Value = moon_vy
Cells(i + 1, 13).Value = moon_px
Cells(i + 1, 14).Value = moon_py

Time = Time + Timestep

i = i + 1

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
If x > 0 And y = 0 Then angle = 90
If x < 0 And y = 0 Then angle = 270
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 = 90 + Atn(y / x) * -1
If y < 0 And x < 0 Then angle = 180 + Atn(x / y)
If y < 0 And x > 0 Then angle = 270 + Atn(y / x) * -1

End Sub

So what this does, step by step:
1. Calculate the moon's angle relative towards the earth using simple trigonometry.
2. Calculate the moon's acceleration for both the X and the Y axis, using trigonometry and the gravitational acceleration equation.
3. Apply acceleration to velocity.
4. Apply velocity to position.
5. Print all variables.
6. Advance time by 1 second.
7. Loop.

This doesn't work.

(At first I even planned to do an RK4 implementation but that was a complete nightmare.)

A snippet of my printed values:
Code (Text):

Time = 0s:
Angle = 0°
Distance = 385000000m
Acceleration X = 0m/s^2
Acceleration Y = -0,002688982324817m/s^2
Velocity X = 1023m/s
Velocity Y = -0,002688982324817m/s
Position X = 1023m
Position Y = 384999999,997311m

Time = 3600s: (after one hour)
Angle = 0,00956571126379358°
Distance = 385000184,85925m
Acceleration X = -0,0000257216115406931m/s^2
Acceleration Y = -0,00268885671887548m/s^2
Velocity X = 1022,95368785554m/s
Velocity Y = -9,6828746017381m/s
Position X = 3683767,39436805m
Position Y = 384982561,007074m

Now I'm not a rocket scientist, but I do suppose that after one hour, at 1km/s, our moon should have moved along the Y axis a bit more than 17.438km.

After around 7 days the moon should be at Y = 0, at 17.5km/h I doubt it will traverse 385000km

Our X position has moved a lot, however, and extrapolating this we'll see the moon move right indefinitely, whilst only going down a relatively small amount. Like this:

(Black arrow what it is currently doing and brown arrow what it should be doing. I also uploaded it as an attachment for future visitors should this image go down.)

I believe that my calculations are wrong. Or at least the vector implementation of said calculations. (I checked the gravitational acceleration at ground level and it was a steady 9.81m/s^2)

I was wondering if anyone could help me get on track with what I'm doing wrong, or what I should be doing instead. I'm also open to the possibility that my entire idea of how orbital physics should work is off, but if anyone could help me out then that would be greatly appreciated.

#### Attached Files:

• ###### norbit2.png
File size:
3.1 KB
Views:
87
Last edited: May 26, 2014
2. May 26, 2014

### Staff: Mentor

Where is the timestep?

Note: this is the first thing that caught my attention, doesn't mean it is the only one, doesn't mean it is the most important one.

3. May 26, 2014

### xzbobzx

I placed it in my acceleration calculation as that seemed the most logical place to put it:

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

Do I need a second timestep at the position?

4. May 26, 2014

### Staff: Mentor

Yes.

$$a = \frac {d^2x} {dt^2}$$

You need to multiply twice by Δt.

5. May 26, 2014

### xzbobzx

Ah okay, that kind of makes sense.

But now I think about it, my timestep is already 1. (I forgot to add this to the main post, doing so now) I don't think that another multiply by 1 is going to solve this.

Is the concept behind what I'm doing (from what you can tell) correct, by the way? I've toyed with the thought that what I'm doing is just an entirely wrong approach in the first place.

6. May 26, 2014

### tms

I assume that the calls to Cells(i, j) return the contents of the spreadsheet cell at row i and column j.

Perhaps you are running up against a limit of the spreadsheet. Spreadsheets usually have a set maximum number of rows and columns, and at one step per second you will probably quickly exceed those limits.

7. May 26, 2014

### xzbobzx

I've tried ramping the Interations up to a million, and even though it took 10 minutes to get it all calculated I could scroll down aaaallll the way to the bottom.

In fact I just checked: http://office.microsoft.com/en-001/...its-HP010342495.aspx?CTT=5&origin=HP005199291

I was exactly 48,576 iterations shy of hitting the limit. One more zero and I would've been a dead man.

8. May 26, 2014

### tms

The last time I tried using Excel the limit was 65536 rows, I think. Even with the increase, there will only be enough to do about half a revolution of the Moon. Changing the time set to a minute probably won't hurt the accuracy enough to matter.

9. May 26, 2014

### MrAnchovy

Angles in Excel are measured in radians so sub Moon_Angle will return the wrong values unless x and y are both greater than zero (which they are for the first 7 days, so this is not your problem).

If you are assuming that the moon's orbit is circular then you should just use a formula for uniform circular motion. If you want to calculate the moon's orbit then you need to use different and more accurate initial conditions.

10. May 26, 2014

### xzbobzx

Now you mention the radians things: that might actually be a problem. I'll have to look into that at school tomorrow.

And for the moon's orbit, meh :D, I don't really care that the moon isn't flying a perfect circular orbit, I just want a body that adheres to gravitational acceleration properly.
In my final simulation I'll probably just hardcode the moon to follow a set path, I just need to figure out where my eventual spacecraft will end up if I send it hurling to the moon at a variable speed. And for that simulation I need the acceleration equations.
...I think.

11. May 26, 2014

### SteamKing

Staff Emeritus
As a general rule, all trig functions (sin, cos, tan) take radian measure for their argument unless specifically stated otherwise in the documentation, whether it is for Excel, another brand of spreadsheet, or most computing languages. About the only devices which accept degree arguments for trig functions as a default are calculators.

12. May 26, 2014

### AlephZero

It "is" a problem, not "might be" a problem! It might not be the only problem, of course.

You don't need all those tests in your moon angle routine. The "atan2" function will sort it all out for you. http://office.microsoft.com/en-gb/excel-help/atan2-HP005208991.aspx [Broken]

Last edited by a moderator: May 6, 2017
13. May 26, 2014

### xzbobzx

That's a very good rule to follow, I had no idea. Thank you!

Hahaha yeah, I've been shot in the foot by trig functions in programming before. And I didn't know about atan2, that's going to be very useful!

Last edited by a moderator: May 6, 2017
14. May 26, 2014

### tms

Note also that you have set up your directions differently from what the trig functions expect. You have set 0 degrees in the positive y direction, increasing clockwise, while the trig functions expect 0 degrees to be in the positive x direction, increasing counter-clockwise.

15. May 26, 2014

### Staff: Mentor

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. Runge-Kutta might indeed be a bit hard to implement in Excel, but the leap-frog method should be as easy as Euler.

For the kind of simulation you want to run, I would follow MrAnchovy's suggestion and just program a circular orbit.

16. May 26, 2014

### tms

The OP is using VBA (a scripting language for Microsoft Office products), which apparently does not have that function.

Last edited by a moderator: May 6, 2017
17. May 26, 2014

### xzbobzx

I'll have to look into the leapfrog integration. I know Euler is really unstable but it was the only solution I had at hand.

As for the circular orbit. Well yes, I could do that and have my moon perfectly rotate around the earth, but I'm not trying to achieve a perfect round orbit. I'm just using the moon as a staging platform to test if my physics work. You might as well read "moon" as "spaceship".

Well dang.
I can just replace 90 with PI*0.5, 180 with PI, 270 with PI*1.5, and then later change them back to angles to display the angle. Shouldn't be a big deal, I think.

18. May 26, 2014

### AlephZero

So why does the web page say "Applies to: Excel 2003"?

Atan2 has been in every version of Basic, and most other programming languages, for about 50 years now.

(But I admit the last version of Excel that I actually used was from Office 97.)

19. May 26, 2014

### xzbobzx

I don't know what to think anymore.

I'll try it out tomorrow and come back with results!

20. May 26, 2014