Orbit by numeric integration (programming question)

Topper

Addon Developer
Addon Developer
Donator
Joined
Mar 28, 2008
Messages
671
Reaction score
32
Points
43
Hello,

I must have a tomatos on my eyes:

when I try to calculate the next position of an Object on it's orbit (I know there are more exact ways that numeric integration, but it's just for fun) , this funktion will not work corectrly (There a no nice orbits):

Code:
//Calling method:
  private void updateGbodyPositions()
  {
    for (int i = 0; i < (gbodies.size()); i++)
        {                
            GBody body1 = gbodies.get(i);
            if (body1.m > centerstar.m) centerstar = body1;            
            for (int i2 = i+1 ; i2 < (gbodies.size()) ; i2++)
            {
                   GBody body2 = gbodies.get(i2);
                   updateV(body1, body2);
          //....
            }
     }
    }

    private void updateV(GBody body1, GBody body2)
    {        
        VECTOR3 er = new VECTOR3(body1.pos.x - body2.pos.x, body1.pos.y - body2.pos.y, body1.pos.z - body2.pos.z);

        float distance = er.length(); 

        er = er.normalize(); //Normalize the vector to value "1"  


        VECTOR3 f1 = er.mul(GRAV_CONST * body1.m * body2.m / distance * distance);  // body.m  = mass of the body
        VECTOR3 f2 = f1.neg();
        
        VECTOR3 a1 = f1.mul(1/body1.m); //dividide f1 by  the mass of body1 
        VECTOR3 a2 = f2.mul(1/body2.m);
        
        body1.v.add(a1); //addition
        body2.v.add(a2);                     
    }
}

And I'm wundering why the following seems to be working good (previous code, factor = gravitational constant). I got this by "try and error".
But i never read to divide somethink by r³, it's always r²...

Code:
private void updateV(GBody star1, GBody star2,  double factor)
    {                   
        float dx = star1.pos.x - star2.pos.x;
        float dy = star1.pos.y - star2.pos.y;
        float dz = star1.pos.z - star2.pos.z;
        
        float distxUp2 = dx * dx;
        float distyUp2 = dy * dy;
        float distzUp2 = dz * dz;     
        
        float distance = (float) Math.sqrt(distxUp2 + distyUp2 + distzUp2);
        float distanceUp3 = distance * distance * distance;
        
        if (distanceUp3 == 0)
        {
            return;
        }

        star2.v.x += factor * dx * star1.m / distanceUp3;
        star2.v.y += factor * dy * star1.m / distanceUp3;
        star2.v.z += factor * dz * star1.m / distanceUp3;        
        star1.v.x += factor * -dx * star2.m / distanceUp3;
        star1.v.y += factor * -dy * star2.m / distanceUp3;
        star1.v.z += factor * -dz * star2.m / distanceUp3;   
    }

So my main question is, why is the first method not working?

Used methodsof my VECTOR3 class:

Code:
package gl_sceneProject.src.sceneProject;
class VECTOR3
{
	float x;
	float y;
	float z;		
	
	VECTOR3 (float x, float y, float z)
	{
		this.x = x;
		this.y = y;
		this.z = z;
	}
	
	VECTOR3 ()
	{
		this.x = 0;
		this.y = 0;
		this.z = 0;
	}
	
	public VECTOR3 neg()
	{
		VECTOR3 v3neg = new VECTOR3();
		v3neg.x = -x;
		v3neg.y = -y;
		v3neg.z = -z;
		return v3neg;
	}

	public VECTOR3 mul(float f)
	{
		VECTOR3 v3mul = new VECTOR3();
		v3mul.x = x * f;
		v3mul.y = y * f;
		v3mul.z = z * f;
		return v3mul;
	}
	
	
	void add(float x, float y, float z)
	{
		this.x += x;
		this.y += y;
		this.z += z;
	}
	
	
	
	void add(VECTOR3 v3)
	{
		this.x += v3.x;
		this.y += v3.y;
		this.z += v3.z;
	}	
	
	
	
	public VECTOR3 normalize()
	{
		float a = (float) Math.sqrt(x*x+y*y+z*z);
		x = x / Math.abs(a);
		y = y / Math.abs(a);
		z = z / Math.abs(a);
		return this;
	}
	
	public VECTOR3 normalize(float length)
	{
		float a = (float) Math.sqrt(x*x+y*y+z*z);
		x = x*length / Math.abs(a);
		y = y*length / Math.abs(a);
		z = z*length / Math.abs(a);
		return this;
	}
	
	public float length()
	{
		return (float) Math.sqrt(x*x+y*y+z*z);
	}
	
	public float distance(VECTOR3 vector3)
	{
		VECTOR3 vdiff = new VECTOR3();
		vdiff.x = Math.abs(vector3.x - x);
		vdiff.y = Math.abs(vector3.y - y);
		vdiff.z = Math.abs(vector3.z - z);
		return vdiff.length();
	}
	
	//....
}
 
Last edited:
How is distance defined in the first example? I mean the variable "distance", not the function "er.distance()"
 
Sorry Urwumpe,

I deleted to much by cutting the code ;-)

float distance = er.length();
 
Sorry Urwumpe,

I deleted to much by cutting the code ;-)

float distance = er.length();

Was that really the right position in your code now... er.length() should ALWAYS be 1.0 now.

EDIT: Never mind, fixed. Then you still divide by distance, when you normalize er, so you calculate with one distance less in the first kind of function.
 
In your code, in the first updateV function, you have:

GRAV_CONST * body1.m * body2.m / distance * distance

which seems like you divide by distance then multiply by distance, hence canceling that bit out. I may be wrong about this, but as far as I know, division and multiplication have similar precedence in C[++] and hence are evaluated left to right, as encountered. C[++] will not realize that what you actually want is

GRAV_CONST * body1.m * body2.m / (distance * distance)

Also, in general, use paranthesis. They are your friends. Especially if you have fractions in your expressions.
 
I know there are more exact ways that numeric integration, but it's just for fun ...

Hi Topper. Unfortunately I cannot answer your question, but I am working on pretty much similar problem now as you are. And I wonder, can you explain what methods are more accurate than numeric integration ?
Maybe if you neglect air drag and other similar perturbation, some analytical method might work, but generally speaking, what can be more accurate than numerical integration ? Thanks !
 
Hi.

I just note that using

VECTOR3 f1 = er.mul(GRAV_CONST * body1.m * body2.m / (distance * distance));

f1 will always be 0, if one of the masses is 0...

In my mind, this is not correct
So this function will also work for zero-masses. In this case, the massiv object will affect the 0-mass-object, but the 0-mass-object will never affect a massiv object.

Code:
    private void updateV2(GBody body1, GBody body2)
    {        
        VECTOR3 v3a = new VECTOR3 (body1.pos);        
        VECTOR3 v3distance = v3a.minus( body2.pos);                
        float distance = v3distance.length();
        float distanceUp3 = distance * distance * distance;        
        if (distanceUp3 == 0) return;
        VECTOR3 vh = v3distance.mul(SystemParameters.GRAV_CONST/distanceUp3);
        body2.v.add(vh.mul(body1.m));
        body1.v.add(vh.mul(-body2.m));
    }
It's the same as solutin "b" in my first post, but using the VECTOR3 class.
It seems, that it is working corect, but i have to explaination why it is corect...
 
Last edited:
well, Einstein needed to develop the general theory of relativity to fix this problem with zero resting mass particles and you just make a dirty hack....
 
I just note that using

VECTOR3 f1 = er.mul(GRAV_CONST * body1.m * body2.m / (distance * distance));

f1 will always be 0, if one of the masses is 0...

In my mind, this is not correct

In Newtonian gravity, it is correct. If something has no mass, there's no gravity that can work on it. That said, I can see why you'd say something like Sputnik has 0 mass when you're otherwise simulating Jupiters.

So this function will also work for zero-masses. In this case, the massiv object will affect the 0-mass-object, but the 0-mass-object will never affect a massiv object.

There are hax to solve this, one of which is you don't count the force caused by the small object on the massive one.

Another is indeed what you propose-

Code:
    private void updateV2(GBody body1, GBody body2)
    {        
        VECTOR3 v3a = new VECTOR3 (body1.pos);        
        VECTOR3 v3distance = v3a.minus( body2.pos);                
        float distance = v3distance.length();
        float distanceUp3 = distance * distance * distance;        
        if (distanceUp3 == 0) return;
        VECTOR3 vh = v3distance.mul(SystemParameters.GRAV_CONST/distanceUp3);
        body2.v.add(vh.mul(body1.m));
        body1.v.add(vh.mul(-body2.m));
    }

Why it works:

F=m*a

So then, what happens if F = G*M*m/(d^2)?

You have (for body m) a = G*M/(d^2)

In other words, the (gravitational) acceleration on body m does NOT contain the number m. Nice.
 
Last edited:
So here is the final function which seems to work perfect :thumbup:

I testet it with original solar system data and it seems to be working.

"factor" has to be used for simulation speed (and accuracy)

Code:
    private void updateV3(GBody body1, GBody body2, float factor)
    {     
    	VECTOR3 d = body1.pos.minus(body2.pos);    	
        float distance = (float) d.length();
        float distanceUp3 = distance * distance * distance;        
        if (distanceUp3 == 0) return;     
        float h = factor * SystemParameters.GRAV_CONST / distanceUp3;      
        body2.v.add(d.mul(h * body1.m));
        body1.v.add(d.mul(h * body2.m));               
    }
 
Last edited:
Now, I want to replace the methods by an RK4 method (Runge-Kutta) and have some questions regarding this because until now, I don't understand the RK4 method by 100%...

So my first question regarding this is, do I have to "Re - Calculate" the forces / accelerations of the attraction between each boddy and each other boddy after eacht step of the RK4 method for the n-body problem? Because if I understand corect, I have a new position and velocity after each step of the RK4 method. If I have another position, the distance will change for sure. And if the distance change, the force and acceleration will also change...?

Hope someone can understand this complicate question...

Until now, I have something like this

Code:
	public static double deriv(int step, double pos, double vel, double t)
	{
		/*
		 * something to calculate
		*/
	}
	
	public static void step (double t, double dt, VECTOR3 pos, VECTOR3 v) 
	{
		double k1_var, k1_vel;
		double k2_var, k2_vel;
		double k3_var, k3_vel;
		double k4_var, k4_vel;
		double[] vel = {v.x,v.y,v.z};
		double[] var = {pos.x,pos.y,pos.z};
		
		for (int i=0; i < var.length; i++) 
		{
			k1_var = vel[i] * dt;
			k1_vel = deriv (i,var[i],vel[i],t)*dt;

                        /*Recalculate force for gravitation???*/

			k2_var =  (vel[i] + 0.5*k1_vel)*dt;
			k2_vel = deriv (i,var[i] + 0.5*k1_var, vel[i] + 0.5*k1_vel, t+0.5*dt)*dt;

                        /*Recalculate force for gravitation???*/

			k3_var =  (vel[i] + 0.5*k2_vel)*dt;
			k3_vel = deriv (i,var[i] + 0.5*k2_var, vel[i] + 0.5*k2_vel, t+0.5*dt)*dt;

                        /*Recalculate force for gravitation???*/

			k4_var =  (vel[i] + k3_vel)*dt;
			k4_vel = deriv (i,var[i] + k3_var, vel[i] + k3_vel, t+dt)*dt;

                        /*Recalculate force for gravitation???*/

			var[i] = var[i] + (k1_var + 2.0*k2_var + 2.0*k3_var + k4_var)/6.0;
			vel[i] = vel[i] + (k1_vel + 2.0*k2_vel + 2.0*k3_vel + k4_vel)/6.0;
		}
	}
 
Last edited:
The code looks about right, provided your "deriv" function does calculate the acceleration correctly at the intermediate steps. Note that you shouldn't need to recalculate the full ephemerides for the involved celestial bodies at each of the 4 RK stages:

The values for the first "deriv" call can be re-used from the last call of the previous step time step.

For the two intermediate steps, it's probably sufficient to do a linear interpolation between the two end points. If that isn't enough, a spherical interpolation should really be sufficient (convert planet positions to spherical polar coordinates, and do a linear interpolation of the spherical coordinates).

So you should really only have to do a single ephemerides evaluation for each full time step.
 
Now, I want to replace the methods by an RK4 method (Runge-Kutta) and have some questions regarding this because until now, I don't understand the RK4 method by 100%...

So my first question regarding this is, do I have to "Re - Calculate" the forces / accelerations of the attraction between each boddy and each other boddy after eacht step of the RK4 method for the n-body problem? Because if I understand corect, I have a new position and velocity after each step of the RK4 method. If I have another position, the distance will change for sure. And if the distance change, the force and acceleration will also change...?

Yes, this is true. RK4 requires you to compute forces along a few points of an 'estimated' trajectory. But those points you can already compute from previous stages of the RK4, so all's fine.

HOWEVER, if you want a 'more accurate than explicit Euler' integration scheme for your program, I would suggest something else- the [ame="http://en.wikipedia.org/wiki/Verlet_integration"]Verlet method[/ame].

The reasons for why the class of methods that include Verlet is better than RK4 are a bit arcane (let's just say those methods capture the physics of your problem better), but in particular the Verlet method is also fairly simple:

let x[k] be the position variables at time instant k; let a[k] be the accelerations at time instant k. Let h be the time step length. (And let v[0] be the initial velocity variables.)

then

x[0] given (initial condition)
x[1] = x[0] + v[0]*h + 0.5*a[0]*h*h (from initial conditions, find the very next position)

and thereafter use

x[k] = 2*x[k-1] - x[k-2] + 0.5*a[k-1]*h*h

To be clear(er hopefully), RK4 'should' be more accurate than Verlet, because RK4 is order 4 and Verlet is order 2. HOWEVER, Verlet should not lead to your planets whirling out into space (which I presume is why you made the change in integration scheme; this has happened to me :P), at least not if your timestep is reasonable. Same thing for RK4 anyway, if your time step is too large the solutions will not make sense. But if what you want is for your solutions to be stable if the underlying system is itself stable, and for your solutions to stay close to the real solutions, Verlet should be enough.

If you feel very adventurous, you could try higher order [ame="http://en.wikipedia.org/wiki/Symplectic_integrator"]symplectic integrators[/ame]. But try RK4, try Verlet, either should be fine for you.
 
Last edited:
Thank you all!

Ok I will try the verlet method first...

But I have a little problem...
The following code seems to be work corect, but only if h=1 (step length) and if I don't multiply "hUp2Mult05" by 0.5...
Is the code corect or is there a mistake?

Code:
    private void updateV_Verlet(GBody body1, GBody body2, double factor) // factor = SystemParameters.SIMSPEED
    {             
        VECTOR3 er = new VECTOR3(body1.pos.x - body2.pos.x, body1.pos.y - body2.pos.y, body1.pos.z - body2.pos.z);
        double distance = er.length(); //distance
        er = er.normalize(); //direction        
        VECTOR3 f2 = er.mul(SystemParameters.GRAV_CONST * body1.m * body2.m * factor / (distance * distance)); //force
        VECTOR3 f1 = f2.neg(); //force other direction (for body2)
        VECTOR3 a1 = f1.mul(1/body1.m); //acceleration for body1
        VECTOR3 a2 = f2.mul(1/body2.m); //acceleration for body2        
        //x[1] = x[0] +         
        //v[0]*h + 0.5*a[0]*h*h (from initial conditions, find the very next position)
        final float h = 1; //Must be 1; h > 1 will not work corectly
        float hUp2Mult05 = h * h; //Normaly, this should be multyplied by 0.5...?
        body1.v = (body1.v.mul(h)).plus(a1.mul(hUp2Mult05));
        body2.v = (body2.v.mul(h)).plus(a2.mul(hUp2Mult05));             
    }

...

//Updating position somewere in the code...

		for (GBody star : gbodies)
			star.pos_new.add(star.v.mul(SystemParameters.SIMSPEED));
 
Last edited:
I made a mistake when listing the Verlet formulas, so that's my bad. They should be

x[0] given (initial condition)
x[1] = x[0] + v[0]*h + 0.5*a[0]*h*h (from initial conditions, find the very next position)

and thereafter use

x[k] = 2*x[k-1] - x[k-2] + a[k-1]*h*h

so no 0.5 factor on acceleration for steps beyond x[1].

About h: yes indeed, you would better use smaller h steps, of course ;) h=1 or smaller should do. In general, if you measure time in, I dunno, Frumions or whatever, then you need to make your h less than 1 Frumion.

The reason for this is that the errors of an integration scheme depend on powers of h. Usually, you'd prefer higher order methods (like RK4, where error is roughly proportional to h^4, or Verlet, where error is roughly proportional to h^2). The higher the order of the integrator, the higher the power of h. But then, obviously, if h > 1 your higher order integrator will do worse than Euler, because h^4 > h^2 > h. If h<1, then you have h^4 < h^2 < h so the higher order integrators are better in this case.

Again, you can measure time in whatever you want (seconds, days, months, I dunno), which also means you measure velocities in (distance/second, or distance/day etc). But whatever unit of measure you use, the size of h should be below 1 in that unit.

Finally, your code isn't Verlet. For one thing, is that a function where you update velocity? Verlet doesn't need to do that. You can get velocities from Verlet too, but during integration itself they are not computed.

Verlet would be this, for body 1 (note a new member, posOld):

aux = body1.pos;
body1.pos = body1.pos.mul(2) - body1.posOld + acc.mul(h*h);
body1.posOld = aux;

Velocity can be estimated from (body1.pos-body1.posOld).mul(1/h)

If you need to compute velocity during the integration phase, then there is velocity Verlet:

x[k+1] = x[k] + v[k]*h + 0.5*a[k]*h*h (<- and this time the 0.5 really is there :P)
v[k+1] = v[k] + 0.5*(a[k]+a[k+1])*h

note that a[k+1] is the acceleration you compute with the system at the x[k+1] position. So you'd compute x[k+1] first then v[k+1].
 
Hello,

sorry for necro posting but I hope it's ok because I want to continue with this project now and after some tests I don't understand RK4 and it works not good...

The code looks about right, provided your "deriv" function does calculate the acceleration correctly at the intermediate steps. ....

So my biggest question at this time is,

Q1. do I have to recalculate the accelerations in K2 (for example) between the position of first object in K2 and second object in K2 and objekt "n" in K2? Or do I have to calculate the acceleration for all the different intermediate steps (K1 to K4) for the first object first in releationship to the "current position" of alle the other objects?

Because at the moment, I have this "Update Function" for the Euler method:

Code:
private void updateEuler()
    {    
        for (GBody star : gbodies)
        {
            star.pos_new.add(star.v.mul(SystemParameters.SIMSPEED));            
        }        
        
        ArrayList<GBody> removeStars = new ArrayList<GBody>();        
        for (int i = 0; i < (gbodies.size()); i++)
        {                
            GBody body1 = gbodies.get(i);
            if (body1.m > centerstar.m) centerstar = body1;            
            for (int i2 = i+1 ; i2 < (gbodies.size()) ; i2++)
            {
                GBody body2 = gbodies.get(i2);                
                if (!body1.equals(body2))
                {
                    updateV1(body1, body2, SystemParameters.SIMSPEED);    
                    removeStars.addAll(collisionDetection(body1, body2));
                }                
            }            
        }
        gbodies.removeAll(removeStars);    
    }

I do it this way because the force between 2 objects is always the same and I only need the negative value of the force vector and this is the more performant way I think...:

Code:
    private void updateV1(GBody body1, GBody body2, double factor )
    {        
        VECTOR3 er = new VECTOR3(body1.pos.x - body2.pos.x, body1.pos.y - body2.pos.y, body1.pos.z - body2.pos.z);
        double distance = er.length();
        er = er.normalize();    
        
        VECTOR3 f2 = er.mul(SystemParameters.GRAV_CONST * body1.m * body2.m * factor / (distance * distance));  //Force calculation has to be once for both objects after, just *-1 (.neg()) is used
        VECTOR3 f1 = f2.neg();     
        
        VECTOR3 a1 = f1.mul(1/body1.m);
        VECTOR3 a2 = f2.mul(1/body2.m);  
        
        body1.v.add(a1);
        body2.v.add(a2);  
               
    }


Q2. Make it sence to use something like "SIMSPEED" or is it the same as dt in the RK4 method?
 
Last edited:
Hello,

here some news regarding my algorithms...

First, i minimized my euler method. Its working good now, but (as we all know) it's not verry accurately.

By calling this method, all vectors would be updated.
I hope it's possible to read the code without comments...

Code:
    private void updateEuler(double dt1)
    {    
        ArrayList<GBody> removeStars = new ArrayList<GBody>();
        for (int i = 0; i < (gbodies.size()); i++)
        {                
            GBody body1 = gbodies.get(i);
            if (body1.m > centerstar.m) centerstar = body1;            
            for (int i2 = i+1 ; i2 < (gbodies.size()) ; i2++)
            {
                GBody body2 = gbodies.get(i2);                
                updateEulerStep(body1, body2, dt1);        
                removeStars.addAll(collisionDetection(body1, body2));                
            }
        }
        gbodies.removeAll(removeStars);    
    }
The "updateEulerStep" method updates the acceleration vector between two objects:
Code:
    private void updateEulerStep(GBody body1, GBody body2, double dt )
    {        
        VECTOR3 er = new VECTOR3(body1.pos.minus(body2.pos));
        er = er.normalize();            
        VECTOR3 f1 = getForce(body1, body2);
        VECTOR3 f2 = f1.neg();          
        VECTOR3 a1 = getAcceleration(f1, body1.m);
        VECTOR3 a2 = getAcceleration(f2, body2.m);        
        body1.v.add(a1.mul(dt));
        body2.v.add(a2.mul(dt));      
        body1.pos_new.add(body1.v.mul(dt));
        body2.pos_new.add(body2.v.mul(dt));
    }
and the getForce and getAcceleration methods are calculating the acceleration:

Code:
    private VECTOR3 getForce(GBody body1, GBody body2)
    {
        VECTOR3 er = body2.pos_new.minus(body1.pos_new);
        double distance = er.length();
        er = er.normalize();    
        return er.mul(SystemParameters.GRAV_CONST * body1.m * body2.m  / (distance * distance));
    }
    

    private VECTOR3 getAcceleration(VECTOR3 f, double m)
    {
        return f.divi(m);
    }
If I try the heun algorithm (or RK4), the result of a two body system is a spiral... With the euler algorithm it's a perfect circle so I don't know whats wrong here???

Code:
    private void updateHeun(double dt1)
    {
        for (GBody body1 : gbodies)
        {
            double dt2 = dt1 / 2;
            body1.a.reset();
            
            for (GBody body2 :gbodies)
                if (!body1.equals(body2))
                    body1.a.add(getAcceleration(getForce(body1, body2), body1.m));
            
            VECTOR3 rBetween = body1.pos_new.plus(body1.v.mul(dt2));
            VECTOR3 vnew = body1.v.plus(body1.a.mul(dt1));
            VECTOR3 rnew = rBetween.plus(vnew.mul(dt2));
            
            body1.pos_new.set(rnew);
            body1.v.set(vnew);
        }    
    }

Euler:
euler.png

Heun:
heun.png

Any hints?
 
Last edited:
Ha I think I solved the "RK4 issue" myself.
Can please someone tell me if it looks good?

First, there is a method "updateRK4".
Because I don't want to render each step, there is h1 (single step length) and dt (total step length per rendering):

Code:
    private void updateRK4(double h1, double dt)
    {
        
        for (double d=0; d<dt; d+=h1)
        {
            double h2 = h1 *0.5;
            for (GBody body1 : gbodies)
            {
                body1.r1 = body1.pos;
                body1.v1 = body1.v;
                body1.a1 = getNewAcceleration(body1.r1,body1);
            }

            for (GBody body1 : gbodies)
                updateK(body1.r1, body1.v1, body1.r2, h2);
            
            for (GBody body1 : gbodies)
                updateK(body1.v1, body1.a1, body1.v2, h2);
    
            for (GBody body1 : gbodies)
                body1.a2 = getNewAcceleration(body1.r2,body1);
            
            //1->3
            for (GBody body1 : gbodies)
                updateK(body1.r1, body1.v2, body1.r3, h2);
            
            for (GBody body1 : gbodies)
                updateK(body1.v1, body1.a2, body1.v3, h2);
            
            for (GBody body1 : gbodies)
                body1.a3 = getNewAcceleration(body1.r3,body1);
            
            //1->4
            for (GBody body1 : gbodies)
                updateK(body1.r1, body1.v3, body1.r4, h1);
            
            for (GBody body1 : gbodies)
                updateK(body1.v1, body1.a3, body1.v4, h1);
            
            for (GBody body1 : gbodies)
                body1.a4 = getNewAcceleration(body1.r4,body1);
            
            for (GBody body1 : gbodies)
            {
                body1.pos.add (rk4mw(body1.v1,body1.v2,body1.v3,body1.v4,h1));    
                body1.v.add   (rk4mw(body1.a1,body1.a2,body1.a3,body1.a4,h1));
            }
        }
        i++;
    }
updateK is used for updating speed and position:

Code:
    private void updateK(VECTOR3 pos0, VECTOR3 v0, VECTOR3 pos1, double dt)
    {
        //pos1.set(pos0.plus(v0.mul(dt))); //Nicer        
        pos1.x = pos0.x + v0.x * dt;  //Faster?
        pos1.y = pos0.y + v0.y * dt;
        pos1.z = pos0.z + v0.z * dt;        
    }
the method "getNewAcceleration" is calculating the acceleration for a given position (r):

Code:
    private VECTOR3 getNewAcceleration(VECTOR3 r, GBody body1)
    {
        VECTOR3 a = new VECTOR3();
        for (GBody body2 :gbodies)
            if (!body1.equals(body2))
                a.add(getAcceleration(getForce(r, body1.m, body2), body1.m));
        return a;
    }
and finaly, rk4mw is calculating the "RK middle value":

Code:
    VECTOR3 rk4mw(VECTOR3 k1, VECTOR3 k2, VECTOR3 k3, VECTOR3 k4, double dt)
    {
        return new VECTOR3
        (
                (k1.x + k2.x*2 + k3.x*2 + k4.x) * dt / 6,
                (k1.y + k2.y*2 + k3.y*2 + k4.y) * dt / 6,
                (k1.z + k2.z*2 + k3.z*2 + k4.z) * dt / 6
        );
    }
 
For those who are interested (and maybe need such an algorithm for Orbiter one day ;-) ):
I did some optimisations and the attached class is the result.
It would be helpful to me if someone could confirm that this is correct...
However, it seems to be ok...

Code:
package gl_sceneProject.src.sceneProject;

import java.util.ArrayList;

public class CalculatorRungeKutta 
{
    ArrayList<GBody> gbodies = null;
    public double update(double h1, int dt, ArrayList<GBody> gbodies)
    {        
        this.gbodies = gbodies;
        double h2 = h1 * 0.5;
        for (int d=0; d<=dt; d++) 
        {
            resetForce();
            
            int k;
            for (k=0; k<3; k++)
            {
                updateForce(k);
                updateAcceleration(k);
                for (GBody body1 : gbodies) 
                {
                    updateK(body1.r[0], body1.v, body1.r, k, (k==3) ? h1:h2);    
                    updateK(body1.v[0], body1.a, body1.v, k, (k==3) ? h1:h2);
                }
            }
            updateForce(k);
            updateAcceleration(k);
            for (GBody body1 : gbodies)
            {
                body1.r[0].add  (rk4mw(body1.v,h1));    
                body1.v[0].add (rk4mw(body1.a,h1));
            }
        }
        return h1*dt;
    }
    
    private void resetForce()
    {
        for (GBody body :gbodies)
            for (int i=0; i<4; i++)
                body.f[i].reset();
    }
    
    private void updateForce(int k)
    {
        for (int i = 0; i < (gbodies.size()); i++)
        {
            GBody body1 = gbodies.get(i);
            for (int i2 = i+1 ; i2 < (gbodies.size()) ; i2++)
            {
                GBody body2 = gbodies.get(i2);
                body1.f[k].add(Physic.getForce(body1.r[k], body1.m, body2));
                body2.f[k].add(body1.f[k].neg());
            }
        }
    }
    
    private void updateAcceleration(int k)
    {
        for (GBody body : gbodies) body.a[k] = Physic.getAcceleration(body.f[k], body.m);
    }
    
    private void updateK(VECTOR3 pos0, VECTOR3[] v0, VECTOR3[] pos, int k, double dt)
    {    
        pos[k+1].x = pos0.x + v0[k].x * dt;
        pos[k+1].y = pos0.y + v0[k].y * dt;
        pos[k+1].z = pos0.z + v0[k].z * dt;        
    }
    
    private VECTOR3 rk4mw(VECTOR3[] k, double dt)
    {
        double dt6 = dt / 6;
        return new VECTOR3
        (
            (k[0].x + k[1].x*2 + k[2].x*2 + k[3].x) * dt6,
            (k[0].y + k[1].y*2 + k[2].y*2 + k[3].y) * dt6,
            (k[0].z + k[1].z*2 + k[2].z*2 + k[3].z) * dt6
        );
    }
}

Code:
package gl_sceneProject.src.sceneProject;

public class Physic 
{
	static VECTOR3 getForce(VECTOR3 virtualPos1, double m1, GBody body2)
	{
		VECTOR3 er = body2.r[0].minus(virtualPos1);
		double distance = er.length();
    	        er = er.normalize();    
    	        er.multhis(SystemParameters.GRAV_CONST * m1 * body2.m  / (Math.pow(distance, 2)));
                return er;
	}
	
	static VECTOR3 getAcceleration(VECTOR3 f, double m)
	{
		if (m!=0) return f.divi(m);
                else return 0;    //??
	}
	static VECTOR3 getForce(VECTOR3 er, double distance, double m1, double m2)
	{
		return er.mul(SystemParameters.GRAV_CONST * m1 * m2 / (Math.pow(distance, 2)));  //2a. f1
	}
	
	static VECTOR3 getForceVector(VECTOR3 pos1, VECTOR3 pos2)
	{	
		return new VECTOR3(pos1.x - pos2.x, pos1.y - pos2.y, pos1.z - pos2.z).normalize();
	}
}
 
Last edited:
Back
Top