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

  • Thread starter xzbobzx
  • Start date
  • #1
10
0

Homework Statement



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.

Homework Equations



[itex]g = -\frac{GM}{r^2}[/itex]

[itex]g =[/itex] gravitational acceleration
[itex]G =[/itex] gravitational constant
[itex]M =[/itex] earth mass
[itex]r =[/itex] orbit radius

The Attempt at a Solution



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

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

'---------

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:
    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:

n6WyM7D.png

(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.

Thank you for your time!
 

Attachments

  • norbit2.png
    norbit2.png
    3.1 KB · Views: 477
Last edited:

Answers and Replies

  • #2
Borek
Mentor
28,702
3,190
moon_px = moon_px + moon_vx

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
10
0
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.

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
Borek
Mentor
28,702
3,190
Do I need a second timestep at the position?

Yes.

[tex]a = \frac {d^2x} {dt^2}[/tex]

You need to multiply twice by Δt.
 
  • #5
10
0
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
tms
644
17
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
10
0
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.

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
tms
644
17
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
pbuk
Science Advisor
Gold Member
2,269
977
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
10
0
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.

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
SteamKing
Staff Emeritus
Science Advisor
Homework Helper
12,796
1,670
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
AlephZero
Science Advisor
Homework Helper
6,994
293
Now you mention the radians things: that might actually be a problem. I'll have to look into that at school tomorrow.

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:
  • #13
10
0
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.

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


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]

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:
  • #14
tms
644
17
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
DrClaude
Mentor
7,601
3,997
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.
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
tms
644
17
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]
The OP is using VBA (a scripting language for Microsoft Office products), which apparently does not have that function.
 
Last edited by a moderator:
  • #17
10
0
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.

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".


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

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
AlephZero
Science Advisor
Homework Helper
6,994
293
The OP is using VBA (a scripting language for Microsoft Office products), which apparently does not have that function.

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
10
0
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.

I don't know what to think anymore.

I'll try it out tomorrow and come back with results!
 
  • #21
pbuk
Science Advisor
Gold Member
2,269
977
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.

You're missing the point. If you want to calculate the actual trajectory of the moon you need to start with the position and velocity of the moon at some point on the moon's actual orbit. As you don't have these you would be better off assuming that the moon has a uniform circular orbit.

As for the trig, it is unneccessary anyway. Note that
$$ F_x = -f cos \theta, F_y = - f sin \theta $$ and
$$ x = r cos \theta, y = r sin \theta; r^2 = x^2 + y^2 $$ so we have
$$ F_x = -\frac{fx}{r}, F_y = -\frac{fy}{r} $$ or
$$ A_x = -\frac{Gmx}{r^3}, A_y = -\frac{Gmy}{r^3} $$

VBA sucks for this (and most things). If you haven't got anything better I'd use javascript like this - which confirms that your code actually comes up with OK results in the first quadrant at least.
 
  • #22
SteamKing
Staff Emeritus
Science Advisor
Homework Helper
12,796
1,670
The OP is using VBA (a scripting language for Microsoft Office products), which apparently does not have that function.

VBA may not have ATAN2 (I haven't checked), but ATAN2 is one of the functions built into Excel itself.

http://office.microsoft.com/en-us/excel-help/math-and-trigonometry-functions-HP005201253.aspx [Broken]
 
Last edited by a moderator:
  • #23
tms
644
17
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.)
Excel itself has the functions ATAN() and ATAN2(), but VBA only has ATN() (note the spelling difference). Don't ask me why; I avoid Office like the plague. Perhaps VBA has it, but I can't find it, and googling finds people saying it doesn't exist.
 
  • #24
2,015
291


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



The calculation is currect. (at least while the moon is in the first quadrant) The acceleration in the Y-direction: 0,002688982324817m/s^2 doesn't vary much in the first hour.

If we pretend it's constant you get Y = (1/2)at^2 = 17.4 km.
 
  • #25
tms
644
17
You're missing the point. If you want to calculate the actual trajectory of the moon you need to start with the position and velocity of the moon at some point on the moon's actual orbit. As you don't have these you would be better off assuming that the moon has a uniform circular orbit.
I think the OP means to use the circular orbit to test his algorithm before using it to calculate the transfer orbit.
 

Related Threads on I'm trying to simulate moon orbit around earth, my moon is broken.

Replies
5
Views
2K
  • Last Post
Replies
7
Views
752
  • Last Post
Replies
7
Views
612
  • Last Post
Replies
5
Views
19K
  • Last Post
Replies
1
Views
9K
  • Last Post
Replies
6
Views
2K
  • Last Post
Replies
4
Views
798
  • Last Post
Replies
1
Views
3K
  • Last Post
Replies
3
Views
3K
  • Last Post
Replies
4
Views
1K
Top