Page 1 of 2
Calculating gravity in realitime ?
Posted: Sun May 10, 2009 2:54 am
by pdwyer
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
Posted: Sun May 10, 2009 5:51 am
by Hades
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
Posted: Sun May 10, 2009 11:06 am
by jack
Paul, have you searched for n-body simulation?
here's a site that has some open-source codes
http://bima.astro.umd.edu/nemo/
Posted: Tue May 26, 2009 2:19 pm
by Lord
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.
Posted: Tue May 26, 2009 2:30 pm
by pdwyer
thanks, need to find an rar decoder now
Posted: Tue May 26, 2009 2:42 pm
by Lord
pdwyer wrote:thanks, need to find an rar decoder now
How about a 'zipped' program?
http://em.q-soft.ch/640VKP.zip
Posted: Tue May 26, 2009 3:29 pm
by pdwyer
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
;=======================================================================
Posted: Tue May 26, 2009 4:51 pm
by Lord
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.
Posted: Wed May 27, 2009 8:48 am
by idle
It seems to be fine.
Posted: Wed May 27, 2009 10:30 am
by Lord
idle wrote:It seems to be fine.
wich one
Posted: Wed May 27, 2009 11:03 am
by idle
Meaning at Pauls post, it'd be fairly rare to get a nice orbit right away.
Posted: Wed May 27, 2009 3:14 pm
by pdwyer
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
Posted: Wed May 27, 2009 4:11 pm
by pdwyer
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
Posted: Wed May 27, 2009 5:59 pm
by Lord
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
Posted: Thu May 28, 2009 1:25 am
by pdwyer
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
