Gravity

[Cervantes]

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

[beard0]

Nice code cervantes. Although not physically correct gravity mathematics, under the initial conditions you set, it gives a good representation. Nice work.

[Cervantes]

What's wrong with the gravity? Are you referring to the lack of a gravitational constant?

[beard0]

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

[zylum]

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]

[Cervantes]

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:

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.

[beard0]

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!

[Cervantes]

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

[zylum]

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 = F*cos(theta)/m and ay = F*sin(theta)/m

here's an updated version of my program... hopefully it's a bit easier to understand

[Catalyst]

the easiest way to approach this would be todo it all in component form thus elminating angles all together

[Cervantes]

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:

[Tony]

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

[Cervantes]

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.

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.

[rizzix]

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

[Cervantes]

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: