Calculating gravity in realitime ?

pdwyer · Post by **pdwyer** » Sun May 10, 2009 2:54 am

I went back to my star / planet orbit code I wrote a few months back but I'm having problems with one little area and I wanted someone with better math than I to give me an idea.

In the code I have a star with a few planets orbiting it, starting the planets anywhere on the screen and giving them a push in any direction will give them some kind of orbit. I can calculate the acceleration due to gravity on each planet at a given point in time, then add the vectors of the planets movement to the gravity accelleration (based on distance and mass) and see where the planet is.

BUT, it's not quite correct.

since I'm working with time slices which are even, as a planet gets close to a star it speeds up, this means that it the location of the planet moves further per slice and this means that a bigger jump in a straight line was made rather lots of little calculations.

Here's where my math is a little poor still. I know that for curves calculus can be used to calc this down to infinity to see where the planet should move if it wasn't making longer straight line hops but I haven't got that far through my math book.

How (not necessarily in code) would I handle this situation? could anyone point me in the right direction of what I need to be reading about to get this right?

Hopefully I've got the question correct to match the problem I think I'm seeing

Hades · Post by **Hades** » Sun May 10, 2009 5:51 am

Hi.

You seem to have the math already right. What you should now do is measuring the change in course (dot product), check if it is above some threshold, and if yes, recompute that last part with smaller timesteps.
At least that was the way I did that, but I'm neither a math guru nor have researched that thoroughly, so there might be a better way...

Have fun,

Hades

jack · Post by **jack** » Sun May 10, 2009 11:06 am

Paul, have you searched for n-body simulation?
here's a site that has some open-source codes http://bima.astro.umd.edu/nemo/

Lord · Post by **Lord** » Tue May 26, 2009 2:19 pm

I did this type of n-body simulation by myself.
Maybe you wan't to have a look at this:

http://em.q-soft.ch/633VKP.rar

And as Hades already stated: it's a question of timesteps.

pdwyer · Post by **pdwyer** » Tue May 26, 2009 2:30 pm

thanks, need to find an rar decoder now

Lord · Post by **Lord** » Tue May 26, 2009 2:42 pm

pdwyer wrote:thanks, need to find an rar decoder now

How about a 'zipped' program?
http://em.q-soft.ch/640VKP.zip

pdwyer · Post by **pdwyer** » Tue May 26, 2009 3:29 pm

Saved me some time! thanks

You get nice clean elipses, I don't, if you go into the options menu and turn trails you can see I get kind of three petaled flower shapes (more or less depending on distance)

The constants at the top can be fiddled with but there's something wrong with my algorithm or the way I'm calculating vectors...

Any tips on where I'm going wrong?

Code: Select all


Structure Star
    LocX.f
    LocY.f
    VectX.f
    VectY.f
    Mass.l 
    Colour.l
    Size.l
EndStructure

; play with these
#StarCount = 7
#ScrWidth  = 1200
#ScrHeight = 900
#StarSize = 4
#DistScale = 600
#StartingVelocity = 500
#RedGiantMass = 500000000
#PlanetMassDiv = 0.5;0

Declare UpdateStarVectors(StarIdx.l)  ;Idx is to skip
Declare.f GravMag(*Star1.star, *Star2.star)
Declare Main()
Declare ReStart()

Global Dim Stars.Star(#StarCount)
Global MaxDist.f = Sqr(Pow(#ScrWidth,2) + Pow(#ScrHeight,2))
Global ShowTrails.l 
Global UseRedGiant.l = #True ;#False


Main()

Procedure ReStart()

    StartDrawing(ImageOutput(0)) 
        Box(0, 0, #ScrWidth, #ScrHeight,0)
    StopDrawing()
    
    
    ReDim Stars.Star(#StarCount)
    
    For I = 1 To #StarCount
       ; Stars(i)\Mass = Random(999) + 1
        Stars(i)\LocX = (Random(#ScrWidth - 2) + 1) * #DistScale
        Stars(i)\LocY = (Random(#ScrHeight - 2) + 1) * #DistScale
        stars(i)\vectx = #StartingVelocity / 2 - Random(#StartingVelocity) 
        stars(i)\vecty = #StartingVelocity / 2 - Random(#StartingVelocity) 
        Stars(i)\Size = #StarSize
        col.l = Random(255) ;the red'er the planet the heavier
        Stars(i)\Colour = RGB(col,Random(255),Random(255));RGB(Random(255),Random(255),Random(255))
        Stars(i)\mass = Pow(col,2) / #PlanetMassDiv
    Next
    
    If UseRedGiant = #True
        Stars(1)\Size = 10
        Stars(1)\mass =  #RedGiantMass
        Stars(1)\Colour = RGB(255,10,10)
        Stars(1)\LocX = #ScrWidth / 2 * #DistScale
        Stars(1)\Locy = #ScrHeight / 2 * #DistScale
        stars(1)\vectx = 0
        stars(1)\vecty = 0
    EndIf

EndProcedure

;=======================================================================

Procedure Main()
    If OpenWindow(0, 0, 0, #ScrWidth, #ScrHeight, "Rnd Star Field", #PB_Window_SystemMenu | #PB_Window_ScreenCentered | #PB_Window_MinimizeGadget); And CreateGadgetList(WindowID(0))
        If CreateImage(0, #ScrWidth, #ScrHeight, 32)
            ImageGadget(0, 0, 0, #ScrWidth, #ScrHeight, ImageID(0))                      ; imagegadget standard
        EndIf
        
        If CreateMenu(0, WindowID(0))
            MenuTitle("File")
                MenuItem( 1, "Start / Restart")
            
            MenuTitle("Options")
                  MenuItem(2, "Show Trails")
                  MenuItem(3, "Red Giant")      
        EndIf
        
        SetMenuItemState(0,3,1)
        frame.q = 0
        
        Repeat 
        
            frame = frame + 1
            
            StartDrawing(ImageOutput(0)) 
            
            ;clear screen
            If ShowTrails = #False
                Box(0, 0, #ScrWidth, #ScrHeight,0)
            EndIf 
    
            For i = 1 To #StarCount
                UpdateStarVectors(i)
                
                If UseRedGiant = #True
                    Stars(1)\LocX = #ScrWidth / 2 * #DistScale
                    Stars(1)\Locy = #ScrHeight / 2 * #DistScale
                    stars(1)\vectx = 0
                    stars(1)\vecty = 0            
                EndIf
            Next
            
            For I = 1 To #StarCount      
    
                If Stars(i)\LocX < 1 Or Stars(i)\LocY < 1 Or Stars(i)\LocX / #DistScale > #ScrWidth Or Stars(i)\LocY / #DistScale > #ScrHeight
                 
                Else
                    If Stars(i)\Size = 1
                        Box(Stars(i)\LocX/#DistScale, Stars(i)\LocY/#DistScale, Stars(i)\Size ,Stars(i)\Size ,Stars(i)\Colour)
                    Else
                        Circle(Stars(i)\LocX/#DistScale, Stars(i)\LocY/#DistScale, Stars(i)\Size / 2 , Stars(i)\Colour) 
                    EndIf
                EndIf
                
                Stars(i)\LocX = Stars(i)\LocX + Stars(i)\VectX
                Stars(i)\LocY = Stars(i)\LocY + Stars(i)\VectY
                            
            Next
            
            StopDrawing()     
            SetGadgetState(0, ImageID(0))
    
            Repeat 
                event.l = WindowEvent()   
                Select event                  
                    Case #PB_Event_Menu           
                        Select EventMenu()  ; To see which menu has been selected   
                            Case 1 ; Restart
                                ReStart()             
                            Case 2 ; Trails
                                If ShowTrails = #False
                                    ShowTrails = #True
                                    SetMenuItemState(0,2,1)
                                Else
                                    ShowTrails = #False
                                    SetMenuItemState(0,2,0)
                                EndIf   
                            Case 3 ; Red giant      
                                If UseRedGiant = #False
                                    UseRedGiant = #True
                                    SetMenuItemState(0,3,1)
                                Else
                                    UseRedGiant = #False
                                    SetMenuItemState(0,3,0)
                                EndIf                                                  
                        EndSelect
                EndSelect               
            Until event = 0 Or event = #PB_Event_CloseWindow 
            If frame % 10 = 0
                Delay(10)
            EndIf
             
        Until event = #PB_Event_CloseWindow                 
    EndIf
EndProcedure

;=======================================================================

Procedure.f GravMag(*Star1.star, *Star2.star)

    Distx.l = (*Star2\LocX - *Star1\LocX)
    Disty.l = (*Star2\LocY - *Star1\LocY)
    
    If Distx < 0: Distx = Distx * -1: EndIf
    If Disty < 0: Disty = Disty * -1: EndIf
    
    Dist.d = Sqr(Pow(Distx,2) + Pow(Disty,2))
    
    ;ProcedureReturn (1 / Pow((Dist),2)) * *Star2\mass
    ProcedureReturn *Star2\mass * 0.006673 / Pow((Dist),2)

EndProcedure

;=======================================================================

Procedure UpdateStarVectors(StarIdx.l)  ;Idx is to skip

    For i = 1 To #StarCount
        If i <> StarIDx 
            ;stop stars shooting off too fast
            If (Stars(i)\locX - Stars(StarIdx)\locX > 1 Or Stars(StarIdx)\locX - Stars(i)\locX > 1) And (Stars(i)\locy - Stars(StarIdx)\locy > 1 Or Stars(StarIdx)\locy - Stars(i)\locy > 1) 
             
        
            ;Get Distance for grav strength
                Magnitude.d = GravMag(@Stars(StarIdx), @Stars(i))
                ;get vector describing diffrence between star 1 and 2 and add
                VectX.f = VectX + ((Stars(i)\locX - Stars(StarIdx)\locX) * Magnitude)
                VectY.f = VectY + ((Stars(i)\locY - Stars(StarIdx)\locY) * Magnitude)

                   
            EndIf 
        EndIf
    Next     

    Stars(StarIdx)\VectX = Stars(StarIdx)\VectX + (VectX * 0.01)
    Stars(StarIdx)\VectY = Stars(StarIdx)\VectY + (VectY * 0.01)
                
EndProcedure

;=======================================================================

Lord · Post by **Lord** » Tue May 26, 2009 4:51 pm

As I understand your code, you are calculating the effect of gravity only in one direction. Let's say you have a sun and a planet. Gravity pulls the planet toward the sun, but the planet is also pulling the sun towards the planet.

So you have to iterate twice through your set of objects. Maybe this will solve your problem.

Another thing:

Code: Select all

    If Distx < 0: Distx = Distx * -1: EndIf
    If Disty < 0: Disty = Disty * -1: EndIf
   
    Dist.d = Sqr(Pow(Distx,2) + Pow(Disty,2))

The square of a digit is always positiv. So there is no need for testing for negativ numbers.

If you are interested, I can give you the source of my program. 1000 lines of code a bit to long for posting and is not good example of good programming.

Post by **idle** » Wed May 27, 2009 8:48 am

It seems to be fine.

Lord · Post by **Lord** » Wed May 27, 2009 10:30 am

idle wrote:It seems to be fine.

wich one

Post by **idle** » Wed May 27, 2009 11:03 am

Meaning at Pauls post, it'd be fairly rare to get a nice orbit right away.

pdwyer · Post by **pdwyer** » Wed May 27, 2009 3:14 pm

The code does go through the stars twice, one loop calling the update proc for each star and then another loop in that proc calling calculating the vector change for each other star's effect (if you turn up the wieght of all the other stars then everything chaotically orbits everthing.

Watching the stars with no trails the movement looks good untill the trails are on you can't see what's happening.

If you look at the orbit of Haleys comet its like in Lords code but in mine the start doesn't do that tight curve around the center star! It's not a distance issue or a mass issue, and it's based on the square of the distance but something is wrong, I need the curve to tighten without the velocity drop I think

need to think about it

pdwyer · Post by **pdwyer** » Wed May 27, 2009 4:11 pm

I think I'm using vectors wrong now that I take a closer look at how the numbers change. It doesn't really use a vector magnitude correctly... I think that may be the problem

Lord · Post by **Lord** » Wed May 27, 2009 5:59 pm

Maybe it helps if I post my calculation routine:

Code: Select all

Procedure Calc()
Protected i.i = 0, j.i=0, dx.d=0, dy.d = 0, dz.d = 0, m2.d = 0, fx.d = 0, fy.d = 0, mm.d = 0

With Obj(i)   While \s <> -1 ; *** Outer loop starts here
    If \s = 1
      NObj(i)\x  = \vx*dt + \ax*0.5*dt*dt + \x; Position and velocity after dt
      NObj(i)\y  = \vy*dt + \ay*0.5*dt*dt + \y
      NObj(i)\vx = \vx    + \ax*dt
      NObj(i)\vy = \vy    + \ay*dt
      fx = 0
      fy = 0
      j  = 0
      While Obj(j)\s <> -1; *** Inner loop starts here
        If \s = 1; object still active
          If i <> j; not the same object
            dx = Obj(j)\x - \x
            dy = Obj(j)\y - \y
            dz = Sqr(dx*dx+dy*dy)
            If dz < (\r + Obj(j)\r); collision of objects
              mm = \m + Obj(j)\m
              NObj(i)\vx = ( \m * \vx + Obj(j)\m * Obj(j)\vx ) / mm; impulse is taken
              NObj(i)\vy = ( \m * \vy + Obj(j)\m * Obj(j)\vy ) / mm
              NObj(j)\vx = NObj(i)\vx; not elastic 
              NObj(j)\vy = NObj(i)\vy
              If (\m >= Obj(j)\m) And Obj(j)\s=1; bigger object 'eats' smaller object
              	Beep_(400,25)
                Objects-1
                StatusBarText(#DataStatusBar, 1, " Active objects: "+Str(Objects))
                NObj(i)\m = mm
                NObj(j)\m = 0
                NObj(i)\r = Sqr(\r*\r + Obj(j)\r*Obj(j)\r); new size
                NObj(j)\r = 0
                NObj(i)\x = (Obj(j)\x * Obj(j)\r + \x * \r) / ( \r + Obj(j)\r ); new position after hit
                NObj(i)\y = (Obj(j)\y * Obj(j)\r + \y * \r) / ( \r + Obj(j)\r )
                NObj(j)\s = 0
                Obj(j)\s  = 0
              EndIf
            Else
              m2 = \m * Obj(j)\m / (dz*dz*dz); sum up effects on object
              fx = fx + m2 * dx
              fy = fy + m2 * dy
            EndIf
          EndIf
        EndIf
        j+1
      Wend
      If \s; still active
        NObj(i)\ax = fx / \m; apply effects
        NObj(i)\ay = fy / \m
      EndIf
    EndIf
    i+1
  Wend
EndWith

Obj.Datas()=NObj.Datas()

EndProcedure ;calc()

I must say, that the underlaying formulas are not developed by myself.
I took them some years ago from an article in 'mc', a german computer magazine.
dt is the chosen time step
Obj() and NObj are structured arrays.
Obj() holds present object datas, NObj() holds the new calculated data.
After all calculations are finished, Obj() becomes NObj().
s is the status of that specific object
x,y are postions
vy,vy are velocities in x- / y-direction
ax,ay are accelerations in x- / y-direction
r is the radius
m is the mass

pdwyer · Post by **pdwyer** » Thu May 28, 2009 1:25 am

Thanks, I'll take a look at these when I get home tonight. I think I know what I'm doing wrong now seeing as the outer part of the ellipse is correct and the inner part not as it speeds up, I'm taking a shortcut with time and velocity. If I take the short cut out I'll probably get the correct ellipse but you won't see any velocity increase for acceleration and deceleration unless I add more complexity to make it correct.

I'm interested to see how yours works as I want acceleration in there too, not just correct ellipses