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

  • Thread starter xzbobzx
  • Start date
  • #26
tms
644
17
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.
 
  • #27
10
0
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 Radian 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
    Radian = 0.0174532925
    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
  • #28
BvU
Science Advisor
Homework Helper
14,401
3,713
Made a few small changes to your code:

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 !:frown:
 
  • #29
BvU
Science Advisor
Homework Helper
14,401
3,713
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 !
 
  • #30
adjacent
Gold Member
1,549
63
Is excel this powerful? I though it's something accountants use for their work. -_-
 
  • #31
10
0
Enjoyed playing with this, even though I just noticed I was a little too late !:frown:

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!
 
  • #32
Borek
Mentor
28,706
3,196
You want to update the vx and vy using the * timestep too !

He did, read the code again.
 
  • #33
BvU
Science Advisor
Homework Helper
14,401
3,713
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...
 
  • #34
BvU
Science Advisor
Homework Helper
14,401
3,713
Is excel this powerful? I though it's something accountants use for their work. -_-
Excel on its own is already extremely powerful. With VBA as programming language with quite some object orientation is powerful squared.

Biggest drawback is that in practice it is the epitome of unstructured programming, so re-use and knowledge transfer are disaster areas.

Having said that, there's an awful lot of science and technology it can be (and is) used for.
 

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
799
  • Last Post
Replies
1
Views
3K
  • Last Post
Replies
3
Views
3K
  • Last Post
Replies
4
Views
1K
Top