perfect orbits by vector math ?

[zillion42]

Hi everyone,

first post... Well, I've done something with the irrlicht engine and without writing too much it's probably best to watch and quote what I've already written..

I made a solar system that simulates newtons and keplers laws of motion using the irrlicht 3d engine... The planets position is calculated according to real data using code byStephen R.Schmitt that I ported from zeno to c/c++:
http://home.att.net/~srschmitt/planetorbits.html

Mean Orbital Elements (J2000)
(Epoche = J2000 = 1.5 January 2000)

Planet a(AU) e i[°] Omega[°] ~omega[°] L[°]
--------------------------------------------------------------------------------
Merkur 0.38709893 0.20563069 7.00487 48.33167 77.45645 252.25084
Venus 0.72333199 0.00677323 3.39471 76.68069 131.53298 181.97973
Erde 1.00000011 0.01671022 0.00005 -11.26064 102.94719 100.46435
Mars 1.52366231 0.09341233 1.85061 49.57854 336.04084 355.45332
Jupiter 5.20336301 0.04839266 1.30530 100.55615 14.75385 34.40438
Saturn 9.53707032 0.05415060 2.48446 113.71504 92.43194 49.94432
Uranus 19.19126393 0.04716771 0.76986 74.22988 170.96424 313.23218
Neptun 30.06896348 0.00858587 1.76917 131.72169 44.97135 304.88003
Pluto 39.48168677 0.24880766 17.14175 110.30347 224.06676 238.92881

The scale of the planets is proportional to each other, the distance is 1/100000.

The simulation of gravity however is greatly simplified and follows the following approach... instead of particles think camera, instead of 1 black hole think 9 planets...

http://www.linuxdevcenter.com/pub/a/linux/2000/09/15/blackhole.html?page=4
The math needed to calculate accurately the motion of objects within a gravity field (in one, two or even three dimensions) is actually quite simple. We'll be using an iterative method to solve what's called the "two body" problem. We simplify the problem even further to assume the smaller particle has no effect on the black hole.

For each moment in time, we update each particle's position based on its current velocity vector. Next, we calculate the square of the distance between our particle and the black hole. We use this as the (possibly scaled) denominator of a gravity constant to find the force of gravity being felt by the particle, thus implementing the inverse-square-law of gravity (and most other forces). The particle's current velocity vector is then updated to include the force of gravity (the vectors are added together). Lastly, we trim a small amount of velocity away from each particle, in order to simulate a slow decay of the orbits.

This simple algorithm produces true orbital dynamics which subscribe to Kepler and Newton's laws, providing the particles don't approach the singularity too closely. As seen in the image above, the orbits are ellipses, as first realized by Johannes Kepler in the early 1600s after studying the orbits of the planets. Newton later came up with his law of universal gravitation, and showed that orbits can also be parabolic or hyperbolic. Both of the latter orbit forms are "open" (the particle escapes), and can appear in our simulation.

The only difference is that I'm checking the gravity pull for 9 planets instead of one black hole (direction and force) and then add all vectors together... vectors with a small length have no effect, that is the same as to say distant planets have a negliable effect...

So to start with the what I try to achieve...
I wanted to make it a little more realistic... So I started with implementing real values not arbitrary game-play tuned values... So if you manually feed the formula below with the distance of earth radius you will get a gforce value of 9.79 something... quite close, and I double checked... So what does this tell me...

void handleG()
{
    pull = core::vector3df(0,0,0);
    u32 i;
    double force = 0;
    for(i=0; i<planets.size(); i++)
    {
        //vector length is strength of gPull
        //direction is determined by vector subtraction
        core::vector3df planPos = planets[i]->getPosition();
        camPos = camera->getPosition();
        core::vector3df cam2plan = planPos-camPos;
        double pdistance = planPos.getDistanceFrom(camPos);
        //distance in meters
        double distanceM = pdistance*scaleF*1000;
        //force in                                        km/s
        //pmass[] is an array with planets mass
        gforce = (6.67428e-11*pMass[i]/pow(distanceM,2))/1000;
        //force in scale factor of sim
        force = gforce/scaleF;
        //set Direction Vector's length to pull
        cam2plan.setLength(force);
        //add the vector for this planet
        //to the total gravity vector.
        pull += cam2plan;
        //printf("%fn",force);
        //printf("%f ,%f ,%f ,n",pull.X,pull.Y,pull.Z);
        if(i==9){
            sunDistanceKM = (pdistance*scaleF);
            sunDistanceAU = (sunDistanceKM/149598000.0);
        }
    }
    //pull -= pull/200;
}

When my camera is on earth surface a 9.81 m/s^2 force pulls it into direction of its center. Let's say my camera is moving tangential to earth surface with 9.81 m/s, should it stay in orbit ?
Of course not, and a stone travelling at 35,316 Km/h tangential to earth surface on sea level would also not...
but somehow thats what this falsely interpreted picture put into my brain...

I got two more videos uploading...

[ame="http://www.youtube.com/watch?v=X2-GdGju5t0"]YouTube - eccentric[/ame]
They show nice orbits, but the values are wrong...
Anyway in a very unscientific way to say it, my thought was that if something is travelling with the same force forward as it's pulled towards the side, it would be a stable orbit... unfortunatley that is not the case... the length(force) of the velocity vector is nearly always longer, except maybe sometimes in perihelion of a very eccentric orbit, than the force of the pull vector, at least according to my observations.

Well to make it shorter, the real values are not working, I have to set my gravity constant to 6.67428e-7 instead of 6.67428e-11 to get good orbits... And I have no clue why...
I'll post the whole source... and continue... I'n not even able to phrase the question properly...
To give it a try... what I do is adding the pull vector to the velocity vector to and the result of that to my position to get my next position. I assumed that if pull vector and velocity vector had same length and were exactly perpendicular to each other I should find myself in a perfect circular orbit, but I was wrong...
Anyone know why ?
Current source is here...
DOWNLOAD
thx

 

[Urwumpe]

Which system do you use for numerically solving the propagation? Runge-Kutta?

 

[zillion42]

the last time I heard the term runge-kutta was when switching dynamics solver models in a 3d application.
To break it down again. There is 2 parts.

1. Planets position is calculated by the probably 'normal way'

1. Find the day number or time since the date of the elements
2. Find the average (or mean) orbital elements of the planet
3. Find the true anomaly, the angle to the planet from perihelion.
4. Find the heliocentric radius, distance to planet from sun
5. Find the heliocentric ecliptic rectangular coordinates of the planet

2. Gravity simulation for camera (<-where I'm stuck)

1. My satellite has a position
2. I have a vector representing my satellites linear motion (direction & speed) pointing in the direction of flight and having a length equal to its speed per frame
3. I have a vector for gravitational pull pointing in the direction of the biggest pull with a length equal to G*planetMass/radius^2... (and like I said calculating 6.67428e-11*330e21/pow(6378.0f,2) yields 9.79something)
4. Each frame I add the velocity vector to the gPull Vector (x1+x2)(y1+y2)(z1+z2) and add the result of that to the satellites position.
5. The result of that becomes my next position, and I start over again...

Surprisingly for me was that if velocity vector and pull vector were perpendicular and of equal length the resulting orbit was far from circular, as I would have expected...

I think that is as clear as I can say it today, I can hardly read anymore...
:tumbleweed:
|
|
\/
no Idea what you mean... the two vectors are the green one in the video (velocity) and the red one (gPull). What is the radius vector and why do I need it ?

 

[Urwumpe]

Surprisingly for me was that if velocity vector and pull vector were perpendicular and of equal length the resulting orbit was far from circular, as I would have expected...

Doh!

Of course not. What was the radius vector then?

 

[zillion42]

It's all fine now

http://irrlicht.sourceforge.net/phpBB2/viewtopic.php?t=32299