Computer Science Canada

Gravity

Author:

Cervantes [ Sat Nov 27, 2004 10:47 am ]

Post subject:

Gravity in 20 lines:

code:

var winID := Window.Open ("position:centre;centre,graphics:200;200,offscreenonly,nobuttonbar,title:Gravity")
var x : array 0 .. 1 of real := init (50, -50)
var y : array 0 .. 1 of real := init (50, -50)
var vx : array 0 .. 1 of real := init (-1.75, 1.75)
var vy : array 0 .. 1 of real := init (0, 0)
var mass : array 0 .. 1 of real := init (1500, 1500)
loop
vx (0) += ((x (1) - x (0)) / Math.Distance (x (0), y (0), x (1), y (1))) * ((mass (0) * mass (1)) / (Math.Distance (x (0), y (0), x (1), y (1)) ** 2)) / mass (0)
vy (0) += ((y (1) - y (0)) / Math.Distance (x (0), y (0), x (1), y (1))) * ((mass (0) * mass (1)) / (Math.Distance (x (0), y (0), x (1), y (1)) ** 2)) / mass (0)
vx (1) += ((x (0) - x (1)) / Math.Distance (x (0), y (0), x (1), y (1))) * ((mass (0) * mass (1)) / (Math.Distance (x (0), y (0), x (1), y (1)) ** 2)) / mass (1)
vy (1) += ((y (0) - y (1)) / Math.Distance (x (0), y (0), x (1), y (1))) * ((mass (0) * mass (1)) / (Math.Distance (x (0), y (0), x (1), y (1)) ** 2)) / mass (1)
cls
for i : 0 .. upper (x)
x (i) += vx (i)
y (i) += vy (i)
drawfilloval (round (x (i) + (maxx / 2)), round (y (i) + (maxy / 2)), 5, 5, black)
end for
View.Update
delay (10)
end loop

If not for having to spend five lines on variable initialization, I could have done gravity between n bodies in 20 lines Confused

Author:	beard0 [ Sat Nov 27, 2004 12:47 pm ]
Post subject:
Nice code cervantes. Although not physically correct gravity mathematics, under the initial conditions you set, it gives a good representation. Nice work.

Author:	Cervantes [ Sat Nov 27, 2004 1:43 pm ]
Post subject:
What's wrong with the gravity? Are you referring to the lack of a gravitational constant?

Author:

beard0 [ Sun Nov 28, 2004 6:18 pm ]

Post subject:

The formula that you use is correct, and the lack of gravitational constant makes sense in that you would also be integrating a pixels/meter constant that would cancel it. However, if you change the initial values such that both objects are at rest, they fly towards each other, pass through each other, and fly off into deep space in opposite directions. For 20 lines your program is excellent. If you were to adapt it to be more than 20 lines, and more correct, some changes that would be necessary to increase the correctness (I say increase because to get it quite accurate, would be impossible) would be:

Decide on the elasticity of collisions between the objects, and implement it into the code.
Calculate acceleration more than once per second, in fact, as close to all the time as possible. This because if one object is flying by the second, very close to it at high speed, it can pick up more speed, and end up well beyond the second before forces and acceleration are re-calculated, which really skews things.

Inspired by your work, I did try to expand on it. For collisions, I decided to have the program end as soon as a collision occurred, so that I wouldn't have to either bounce the objects off of each other, or figure out the spinning path that they would take if they were joined. And to increase the validity of the accelerations, I calculated them, moved the object one unit in the right direction, and added the proportional amount to a time variable in a loop, exiting once the time reached 1 second. Here's my program, if you're interested. I don't believe it to be completely accurate, simply as accurate as I could make it without spending too much time on it. Note: This is an expansion on Cervantes' work.

code:

var winID := Window.Open ("position:centre;centre,graphics:max;max,offscreenonly,nobuttonbar,title:Gravity")
var x : array 0 .. 1 of real := init (50, -50)
var y : array 0 .. 1 of real := init (50, -50)
var vx : array 0 .. 1 of real := init (-1.75, 1.75)
var vy : array 0 .. 1 of real := init (1.75, -1.75)
var mass : array 0 .. 1 of real := init (1500, 1500)
var r : array 0 .. 1 of int
if mass (0) < mass (1) then
r (0) := 5
r (1) := round (5 * mass (1) / mass (0))
else
r (1) := 5
r (0) := round (5 * mass (0) / mass (1))
end if
proc accelx
vx (0) += ((x (1) - x (0)) / Math.Distance (x (0), y (0), x (1), y (1))) * ((mass (0) * mass (1)) / (Math.Distance (x (0), y (0), x (1), y (1)) ** 2)) / mass (0)
vx (1) += ((x (0) - x (1)) / Math.Distance (x (0), y (0), x (1), y (1))) * ((mass (0) * mass (1)) / (Math.Distance (x (0), y (0), x (1), y (1)) ** 2)) / mass (1)
if Math.Distance (x (0), y (0), x (1), y (1)) <= r (0) + r (1) then
vx (0) := 0
vx (1) := 0
end if
end accelx
proc accely
vy (0) += ((y (1) - y (0)) / Math.Distance (x (0), y (0), x (1), y (1))) * ((mass (0) * mass (1)) / (Math.Distance (x (0), y (0), x (1), y (1)) ** 2)) / mass (0)
vy (1) += ((y (0) - y (1)) / Math.Distance (x (0), y (0), x (1), y (1))) * ((mass (0) * mass (1)) / (Math.Distance (x (0), y (0), x (1), y (1)) ** 2)) / mass (1)
if Math.Distance (x (0), y (0), x (1), y (1)) <= r (0) + r (1) then
vy (0) := 0
vy (1) := 0
end if
end accely
proc newpos
accelx
accely
for i : 0 .. upper (x)
if abs (vx (i)) <= 1 then
x (i) += vx (i)
else
var t := 0.0
loop
x (i) += sign (vx (i))
t += abs (1 / vx (i))
exit when t >= 1
accelx
exit when vx (i) = 0
end loop
end if
if abs (vy (i)) <= 1 then
y (i) += vy (i)
else
var t := 0.0
loop
y (i) += sign (vy (i))
t += abs (1 / vy (i))
exit when t >= 1
accely
exit when vy (i) = 0
end loop
end if
end for
end newpos
colourback (black)
loop
cls
newpos
for i : 0 .. upper (x)
drawfilloval (round (x (i) + (maxx / 2)), round (y (i) + (maxy / 2)), r (i), r (i), 4 + 5 * i)
end for
View.Update
exit when Math.Distance (x (0), y (0), x (1), y (1)) <= r (0) + r (1)
Time.DelaySinceLast (50)
end loop

Author:	zylum [ Sun Nov 28, 2004 6:52 pm ]
Post subject:
your program doesnt seem too smooth... the orbits look sort of rectangular. idk maybe its just me. i know its far from 20 lines but here's my gravity proggy [EDIT] update version is in a following post [/EDIT]

Author:

Cervantes [ Sun Nov 28, 2004 8:37 pm ]

Post subject:

I'd have to agree with you, Zylum. Those orbits do look kinda rectangular. I think it might have something to do with your loops for calculating the positions (x and y).

code:

loop
x (i) += sign (vx (i))
t += abs (1 / vx (i))
exit when t >= 1
accelx
exit when vx (i) = 0
end loop

that, and also the other one for y.
See, you're updating x a bunch of times in that loop, but using the same y value each time you calculate your new x velocity. Then, when you finally exit that loop and go into the y loop, you are calculating the new y velocity over and over based on one x value. The x and y velocities aren't changing at the same time, if you follow me.

nice stuff zylum. I don't know if you've posted this anywhere else (I think I've seen something like this before), but, nonetheless, +25 BITS.

I don't quite understand your math, though. Where you change the accelerations:

code:

BODY (i)._AX += ((_G * BODY (i)._MASS * BODY (j)._MASS) / ((BODY (i)._X - BODY (j)._X) ** 2 + (BODY (i)._Y - BODY (j)._Y) ** 2))
* (cosd (GETANGLE (BODY (j)._X, BODY (j)._Y, BODY (i)._X, BODY (i)._Y))) / BODY (i)._MASS
BODY (i)._AY += ((_G * BODY (i)._MASS * BODY (j)._MASS) / ((BODY (i)._X - BODY (j)._X) ** 2 + (BODY (i)._Y - BODY (j)._Y) ** 2))
* (sind (GETANGLE (BODY (j)._X, BODY (j)._Y, BODY (i)._X, BODY (i)._Y))) / BODY (i)._MASS

Well, the best I can get out of that is this:
Posted Image, might have been reduced in size. Click Image to view fullscreen.

that's for X acceleration.
Why are you dividing by the mass of body 1 when, in the beginning, you multiplied by the mass of body 1. I just took all occurances of body 1's mass (when I say body 1, I mean i, not j) and it worked just the same.

Another thing, shouldn't you be dividing by the square of the distances? not just the xdiff squared + the y diff squared.

Author:	beard0 [ Sun Nov 28, 2004 10:16 pm ]
Post subject:
A good point with the calculating x and y acceleration seperately. In terms of formula, Cervantes has the right one, though I find it funny that you suggest zylum cancel the mass1's, as your formulae actually have mass1*otherstuff/mass1 in them. Zylum, you've done a nice job of user freindliness, though your code has the same problems I outlined in Cervantes. The truth is, the accelerations should be calculated as often as possible. As /\t approaches 0, the physics approaches accuracy . Hmm writing that down makes me think... I wonder if I can use calculus techniques.... allow me to go away and ponder, I shall return with what I come up with!

Author:

Cervantes [ Sun Nov 28, 2004 10:31 pm ]

Post subject:

beard0 wrote:

The truth is, the accelerations should be calculated as often as possible. As /\t approaches 0, the physics approaches accuracy

How inaccurate are our calculations, then?

beard0 wrote:

though I find it funny that you suggest zylum cancel the mass1's, as your formulae actually have mass1*otherstuff/mass1 in them

True.

The reason for that is I just crunched two variables, force and distance, into the four calculations. See below for the more normal version of my code.

It appears as though our calculations agree, Zylum:
Taking your program as a basis, try substituting the following code for your UPDATEBODIES procedure:

code:

proc UPDATEBODIES
for a : 1 .. upper (BODY)
for b : a + 1 .. upper (BODY)
var distance := Math.Distance (BODY (a)._X, BODY (a)._Y, BODY (b)._X, BODY (b)._Y)
var force := _G * (BODY (a)._MASS * BODY (b)._MASS) / (distance * distance)

BODY (a)._VX += ((BODY (b)._X - BODY (a)._X) / distance) * force / BODY (a)._MASS
BODY (a)._VY += ((BODY (b)._Y - BODY (a)._Y) / distance) * force / BODY (a)._MASS

BODY (b)._VX += ((BODY (a)._X - BODY (b)._X) / distance) * force / BODY (b)._MASS
BODY (b)._VY += ((BODY (a)._Y - BODY (b)._Y) / distance) * force / BODY (b)._MASS
end for
end for
for I : 1 .. upper (BODY)
BODY (I)._X += BODY (I)._VX
BODY (I)._Y += BODY (I)._VY
end for
end UPDATEBODIES

and import Math at the top of the class. Or if you don't have 4.0.5, make your own Distance function.

I find that they work identically, though I haven't watched them for all that long So cool

Author:	zylum [ Mon Nov 29, 2004 1:27 am ]
Post subject:
although they may seem similar, they are completetly different. your method would be fine for 2 bodies but since there are more than that, you will first need to calculate the net acceleration before you update the velocity. also i think you need to trig to find the x and y accelerations from the net force... hence the / Mass that beardo pointed out. since F = ma, then a = F/m or in component form, ax = Fcos(theta)/m and ay = Fsin(theta)/m here's an updated version of my program... hopefully it's a bit easier to understand

Author:	Catalyst [ Mon Nov 29, 2004 2:36 am ]
Post subject:
the easiest way to approach this would be todo it all in component form thus elminating angles all together

Author:	Cervantes [ Mon Nov 29, 2004 4:00 pm ]
Post subject:
zylum wrote: your method would be fine for 2 bodies but since there are more than that, you will first need to calculate the net acceleration before you update the velocity. it doesn't? It seemed to act exactly the same to me As for attaction between more than two bodies:

Author:	Tony [ Mon Nov 29, 2004 8:09 pm ]
Post subject:
Cervantes, I don't know if you're aware, but there's a problem with your model.. mainly the trend of the system moving upwards and the fact that one of the moons manages to escape

Author:

Cervantes [ Mon Nov 29, 2004 9:31 pm ]

Post subject:

I believe the general trend of it moving upward is not a problem: it is the entire system (including the sun) moving upwards, not just the planets. Try adding these lines just before the View.Update:

code:

locate (1, 1)
put "SunY : ", cbody (0).y
locate (2, 1)
put "SunYP: ", cbody (0).y * PperD

As you can see, the planets do have an effect on the sun, albeit it is small.
The way I figure, it is because when the program begins, all the planets are lined up on one side. So, the combined mass of all of them begins to pull the sun towards them. As they move around the sun, however, they pull it in different directions. They will never fully negate the beginning tug, however. Bah! I don't quite know how to describe it. Perhaps someone will know what I am talking about and be able to put it into words. Perhaps not.

In sss.t, there are two lines of commented code that look like this:

code:

cbody (0).x := 0 %maintain sun in centre of screen
cbody (0).y := 0

In the attachment, they are commented out. If you uncomment them, the sun stays in the middle and the planets maintain their orbits, for the most part. Saturn and Jupiter do, quite nicely. Earth's orbit is actually quite eccentric. It's aphelion (furtherst point from centre in orbit) and it's perihilion (closest point to centre in orbit) move around. First, it's aphelion may be in the north (at the top of its orbit, in relation to the screen) and its perihelion in the south, then, after time, they gradually reverse themselves, such that the aphelion is in the south and the perihelion is in the north. This, in actuallity, is called the precession of the equinoxes, and it seems to be demonstrated in my model quite well. Smile

There is no denying that Ganymede leaves orbit around Jupiter. One time, even Luna left Earth's orbit. Though when Luna left orbit around Earth, it began to orbit the Sun in a very eccentric orbit. It showed the precession of the equinoxes even more clearly than did Earth!

Wow, and after watching it for even longer, I noticed that Venus' orbit changed from quite circular to quite eccentric, which, once again, actually happens in our solar system.

Aah, and after even longer, Luna was finally ejected from the solar system. I think it must have come to close to Mars or Earth or perhaps both at one time.

This picture shows two instances of a precession of the equinoxes: one is Venus' and the other is Luna's.

Author:	rizzix [ Mon Nov 29, 2004 10:24 pm ]
Post subject:
Cervantes wrote: zylum wrote: your method would be fine for 2 bodies but since there are more than that, you will first need to calculate the net acceleration before you update the velocity. it doesn't? It seemed to act exactly the same to me As for attaction between more than two bodies: I have to agree with zylum.. that will work only for a 2 body system.. having more than 2.. well lets say its just far tooo complex. you would need to find the Net force acting on each planet.. and based on that netforce you would describe that planets elliptical motion

Author:	Cervantes [ Wed Dec 01, 2004 6:49 pm ]
Post subject:
Instead of calculating the net acceleration of each body, I just skipped right to adding to the velocity. And, after the two for loops are done, the net velocity is determined. I admit I don't know that much about this stuff, but perhaps this site does? That's where I got this model. I was playing around with different values for bodies and I always found my model to be more stable than Zylum's. In the following attachments, a three body system is modelled using both my approach and zylum's approach. Given that there are bodies of equal mass at an equal distance from it on opposite sides, should the red body in the centre not remain still? Here's the comparison:

Author:	rizzix [ Wed Dec 01, 2004 8:23 pm ]
Post subject:
whether the centre body remains still depends on the next force acting on it.. which in turn depends on the masses sorrounding it.. but i'm not sure what zylum has done there... i was just reading up your arguments.. and took side of the one that imo is more correct.

Author:	Tony [ Wed Dec 01, 2004 10:35 pm ]
Post subject:
the body in the middle should remain at rest.. though in the code, there's an error that occurs on the 3rd pass in 10^-8 of a unit. That difference in the net force quickly escalates.

Author:	Cervantes [ Thu Dec 02, 2004 4:20 pm ]
Post subject:
So there is. Perhaps this is the kind of problem that beard0 is referring to. It is interesting to note that zylum's approach results in major problems by only the first orbit.

Author:	zylum [ Thu Dec 02, 2004 10:44 pm ]
Post subject:
don't mind me, i was just being a moron my problem was that after i calculated each new acceleration per body, i updated the velocity and position before i finished calculating the other bodies... here's the updated version which does all the calculations FIRST and then updates the positions. now it runs almost the same as yours.. it loses control after about five and a half. i wonder what's still causing the slight differences and which is more accurate?

Author:	rizzix [ Fri Dec 03, 2004 1:28 am ]
Post subject:
pfft i finally ran ur code(s)... heh 3 body system eh? interesting...

Author:	MihaiG [ Wed Dec 22, 2004 3:12 pm ]
Post subject:
it don't work on mine i have 4.0.3x math.distance don't work... teh link for 4.0.5 update doesnt work can some1 post it here... and is there any way around math.distance

Author:	Cervantes [ Wed Dec 22, 2004 7:48 pm ]
Post subject:
sigh I already taught you about how to do it. I think the post got deleted, though, so I'll give you another link. Beard, any updates?

Author:

Cervantes [ Mon Jan 31, 2005 8:45 am ]

Post subject:

In talking on IRC, I got back to thinking about this. A little while later, I had:

Turing:

var x : array 0 .. 3 of real := init (-100, 100, 100, -100)
var y : array 0 .. 3 of real := init (-100, -100, 100, 100)
var vx : array 0 .. 3 of real := init (2.5, 0, -2.5, 0)
var vy : array 0 .. 3 of real := init (0, 2.5, 0, -2.5)
var mass : array 0 .. 3 of real := init (1500, 1500, 1500, 1500)
fcn calculate (n, i : int) : boolean
vx (n) += ((x (i) - x (n)) / Math.Distance (x (n), y (n), x (i), y (i))) * ((mass (n) * mass (i)) / (Math.Distance (x (n), y (n), x (i), y (i)) ** 2)) / mass (n)
vy (n) += ((y (i) - y (n)) / Math.Distance (x (n), y (n), x (i), y (i))) * ((mass (n) * mass (i)) / (Math.Distance (x (n), y (n), x (i), y (i)) ** 2)) / mass (n)
result false
end calculate
loop
exit when calculate (0, 1) or calculate (1, 0) or calculate (0, 2) or calculate (2, 0) or calculate (0, 3) or calculate (3, 0) or calculate (1, 2) or calculate (2, 1) or calculate (1, 3) or
calculate (3, 1) or calculate (2, 3) or calculate (3, 2)
cls
for i : 0 .. upper (x)
x (i) += vx (i)
y (i) += vy (i)
drawfilloval (round (x (i) + (maxx / 2)), round (y (i) + (maxy / 2)), 5, 5, black)
end for
delay (10)
end loop

Gravity between n bodies in twenty lines, thanks to the false function trick!
Mind you, things go berserk after a few orbits.
Bead0, if you're still alive, have you got any further on the program?

You can syndicate this boards posts using the file backend.php or view the topic map using sitemap.php.

Terms of Use | Privacy Policy