----------------------------------- Cervantes Sat Nov 27, 2004 10:47 am ----------------------------------- Gravity in 20 lines: 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 Sat Nov 27, 2004 12:47 pm ----------------------------------- Nice code cervantes. Although not physically correct gravity mathematics, under the initial conditions you set, it gives a good representation. Nice work. :D ----------------------------------- Cervantes Sat Nov 27, 2004 1:43 pm ----------------------------------- What's wrong with the gravity? Are you referring to the lack of a gravitational constant? ----------------------------------- beard0 Sun Nov 28, 2004 6:18 pm ----------------------------------- 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. 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)) = 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: 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: http://www.nhswork.000k2.com/GravityMath.jpg 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 Sun Nov 28, 2004 10:16 pm ----------------------------------- 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 :D. 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 Sun Nov 28, 2004 10:31 pm ----------------------------------- 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? though I find it funny that you suggest zylum cancel the mass1's, as your formulae actually have mass1*otherstuff/mass1 in them True. :lol: 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: 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 :coolest: ----------------------------------- zylum Mon Nov 29, 2004 1:27 am ----------------------------------- 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 Mon Nov 29, 2004 2:36 am ----------------------------------- the easiest way to approach this would be todo it all in component form thus elminating angles all together ----------------------------------- Cervantes Mon Nov 29, 2004 4:00 pm ----------------------------------- 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 :think: As for attaction between more than two bodies: ----------------------------------- Tony Mon Nov 29, 2004 8:09 pm ----------------------------------- Cervantes, I don't know if you're aware, but there's a problem with your model.. http://www.compsci.ca/bbs/files/sss_screen.jpg mainly the trend of the system moving upwards and the fact that one of the moons manages to escape :? ----------------------------------- Cervantes Mon Nov 29, 2004 9:31 pm ----------------------------------- 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: 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: 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. [url=http://www.nhswork.000k2.com/SSS2.jpg]This picture shows two instances of a precession of the equinoxes: one is Venus' and the other is Luna's. ----------------------------------- rizzix Mon Nov 29, 2004 10:24 pm ----------------------------------- 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 :think: 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 Wed Dec 01, 2004 6:49 pm ----------------------------------- 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 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: ----------------------------------- rizzix Wed Dec 01, 2004 8:23 pm ----------------------------------- 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. ----------------------------------- Tony Wed Dec 01, 2004 10:35 pm ----------------------------------- 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. :? http://www.compsci.ca/bbs/files/sss_screen_2.jpg ----------------------------------- Cervantes Thu Dec 02, 2004 4:20 pm ----------------------------------- 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. ----------------------------------- zylum Thu Dec 02, 2004 10:44 pm ----------------------------------- don't mind me, i was just being a moron :lol: 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? ----------------------------------- rizzix Fri Dec 03, 2004 1:28 am ----------------------------------- pfft i finally ran ur code(s)... heh 3 body system eh? interesting... ----------------------------------- MihaiG Wed Dec 22, 2004 3:12 pm ----------------------------------- 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 ----------------------------------- Cervantes Wed Dec 22, 2004 7:48 pm ----------------------------------- *sigh* I already taught you about how to do it. I think the post got deleted, though, so I'll give you another [url=http://www.compsci.ca/v2/viewtopic.php?p=40459]link. Beard, any updates? ----------------------------------- Cervantes Mon Jan 31, 2005 8:45 am ----------------------------------- In talking on IRC, I got back to thinking about this. A little while later, I had: 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?