Orbits

Game Mechanics
#1

Hello, I’m kinda experimenting a little with gravitational forces and orbits in 2D. I have some objects (asteroids or something) being influenced by n static bodies (planet, stars). 2 questions:

  1. What integrator should I use? Euler sucks, it requires too small time steps to be useful, and performance is fairly important. I would really like a code example for gravitation specifically. I’ve found numerous articles on through Google on integration, but I didn’t really understand the actual implementations… T___T
  2. How do I calculate the gravitational acceleration? With the standard G*M/R^2, I just get an acceleration, not an 2D acceleration vector. I currently get the angle between the 2 bodies using atan2(dy, dx) and then get the vector using cos and sin, but it is really slow. I suppose I could compute a direction vector by normalizing (dx, dy), but that would require a sqrt and multiple divides… Is there any faster way? It would be awesome if there was an equation that actually produced a vector in the first place.
    I may be using some words wrong here, so please don’t flame me. xD
#2

  1. Can’t help much, don’t know what integration or Euler is :S

  2. Faster sin/cos using lookup table: http://www.java-gaming.org/index.php/topic,24191.0.html
    Faster atan2 also using lookup table: http://www.java-gaming.org/index.php/topic,24193.0.html

#3

The acceleration from gravity is the length of the acceleration vector, to get the vector, just subtract the position of the object you are being accelerated towards from the object you are accelerating. Then normalise that vector, that gives you a unit vector describing the direction of acceleration, then you can multiply it by the acceleration from gravity and you get the acceleration vector.

then you can just scale a copy of that by time to get the change in velocity for each frame, then apply that to the current velocity vector for the object, then multiply that by time to get the change in position for the frame. Apply to the transform and job done :slight_smile:

Add in angular acceleration + velocity etc and all is nice :).

It’s what I’m doing in darkvoid which is the full 3D 6dof :).

Endolf

#4

Check this thread:
http://www.java-gaming.org/index.php/topic,24085.0.html

There is code for solution with component calculations if you don’t want use vectors. Its quite simple and efficient too. Some cleaning is needed for more general use but you should see how to use it.


public final static long staticGravityForce = 1000; //just test what is good
public final static long staticX = 400;
public final static long staticY = 400;
public final static long staticRadius = 50;
public void gravity(double x, double y, int radius, double dx,double dy) {
      
    double ddx = (x - staticX);
    double ddy = (y - staticY);
    double distanceSquared = ddx * ddx + ddy * ddy;
    double distance = Math.sqrt(distanceSquared);
 
    if ((distance) > (staticRadius + radius)) {
         
        ddx /= distance;
        ddy /= distance;
        dx += (ddx * staticGravityForce / distanceSquared);
        dy += (ddy * staticGravityForce / distanceSquared);             
    } else {
        //COLLAPSE
        //JUST make some cool effects here
    }
    }
}
#5

[quote]1. What integrator should I use? Euler sucks, it requires too small time steps to be useful
[/quote]
I wouldn’t be too quick to dismiss Euler integration. I’d expect it to be perfectly adequate for games (if that’s what you’re doing), and certainly for prototyping. Make the time steps small if you have to, but get everything working before trying to optimize.

But if you’re set on using higher-order methods, try the chapter on “Integration of Ordinary Differential Equations” in Numerical Recipes. Runge-Kutta methods are usually a good starting point.

There presumably are methods that are specific to gravitation, but I’m not familiar with any. (You could maybe try tweaking the kinetic energy of the asteroids so that kinetic plus potential is conserved during each step.)

[quote]2. I suppose I could compute a direction vector by normalizing (dx, dy), but that would require a sqrt and multiple divides…
[/quote]
Normalizing the vector (sqrt) makes more sense than dealing with angles (sin and cos and atan2). Also, knowing the (non-squared) distance of the asteroid from the planet could turn out to be useful in other parts of your simulation.

[quote]I have some objects (asteroids or something) being influenced by n static bodies (planet, stars).
[/quote]
Do the asteroids influence each other gravitationally? If not then the gravitational field of the planets is static, and you can pre-compute it. That is, compute the force vector (or potential value) at each pixel on the screen (or over some other grid) and just look-up the value for the asteroid’s current position (interpolating between grid points if you like).

Simon

#6

I’m using the same solution as Pitbuller. The only thing I have to do is getting rid of atan2. Sqrt should at least be faster. xd
Euler integration is fine with circular orbits, but the whole orbit starts to rotate when the orbit is elliptic. Look at the first test in this article: http://codeflow.org/entries/2010/aug/28/integration-by-example-euler-vs-verlet-vs-runge-kutta/ (run it by holding your mouse over it).

The main goal is to get stable orbits with good performance. Asteroids do NOT influence each other. I may however have multiple stars/planets affecting the body. From a performance perspective Euler integration is perfect, but the problem with high acceleration needs a solution. I’m going to experiment with using a smaller time step when the gravitational force gets higher. That would only increase the computation cost for a fraction of the time. Will report back! :wink:

EDIT: Adapting the time step to the force was a bad idea. The orbit gets unstable…
Does anybody have an implementation of RK4 or the means to teach me how to implement it myself? I would like to test it. xd

#7

It might be worth giving the semi-implicit Euler method a try. It conserves (kinetic plus potential) energy, which might help with the stability of your orbits.
Simon

#8

Here you go. I’ve included a simplified version of one calculation I was using it for as a test case.



/** 
 * Represents a point in R^4. Doubles are used. Instances of this class are immutable.
 * @author Peter Taylor
 */
public class Vec4
{
    private final double w, x, y, z;
    
    /** Constructs a point in R^4 from its projections onto the four axes. */
    public Vec4(double _w, double _x, double _y, double _z) { w = _w; x = _x; y = _y; z = _z; }
    
    /** Vector addition. Adds components in each axis. */
    public Vec4 add(Vec4 q) { return new Vec4(w + q.w, x + q.x, y + q.y, z + q.z); }
    /** Scalar multiplication. */
    public Vec4 scale(double scalar) { return new Vec4(w * scalar, x * scalar, y * scalar, z * scalar); }
    /** Combines scalar multiplication and addition. */
    public Vec4 addScaled(Vec4 q, double scalar) {
        return new Vec4(w + q.w * scalar, x + q.x * scalar, y + q.y * scalar, z + q.z * scalar);
    }
    
    /** Gets projection onto w axis; this is the first component. */
    public double getW() { return w; }
    /** Gets projection onto x axis; this is the second component. */
    public double getX() { return x; }
    /** Gets projection onto y axis; this is the third component. */
    public double getY() { return y; }
    /** Gets projection onto z axis; this is the fourth component. */
    public double getZ() { return z; }
    
    /** Output as "(w, x, y, z)". */
    public String toString() { return "(" + w + ", " + x + ", " + y + ", " + z + ")"; }

    /** Function from (time, Vec4) to Vec4 */
    public static interface Fn {
        Vec4 eval(double t, Vec4 arg);
    }

    /** The most familiar of the fourth order Runge-Kutta methods for first order differential equations. */
    public static class RK4 {
        private final Fn fn;

        public RK4(Fn _fn) { fn = _fn; }

        /** Takes a step in the integration. This is defined by <pre>
         * k1 = h f(t, X)
         * k2 = h f(t + h/2, X + k1/2)
         * k3 = h f(t + h/2, X + k2/2)
         * k4 = h f(t + h, X + k3)
         * X' = X + (k1 + 2k2 + 2k3 + k4)/6</pre>
         * @param X The current position in R^4
         * @param h The step size */
        private Vec4 takeStep(Vec4 X, double t, double h) {
            double h2 = h/2;
            Vec4 k1 = fn.eval(t, X);
            Vec4 k2 = fn.eval(t + h2, X.addScaled(k1, h2));
            Vec4 k3 = fn.eval(t + h2, X.addScaled(k2, h2));
            Vec4 k4 = fn.eval(t + h, X.add(k3));
            Vec4 sum = k2.add(k3);
            sum = k1.addScaled(sum, 2).add(k4);
            return X.addScaled(sum, h / 6);
        }

        /** Solves the equation from t1 to t2 in steps of delta.
         * t2 must be at least as great as t1.
         * @param X The value of X(t1)
         * @param t1 The start value of t
         * @param t2 The end value of t
         * @param delta The step size
         * @return The computed value of X(t2), using dX/dt = fn(X). */
         public Vec4 integrate(Vec4 X, double t1, double t2, double delta) {
             double range = t2 - t1;
             if (range < 0) throw new IllegalArgumentException(t2 + " is smaller than " + t1 + " - negative range");
             // Do as much as can be done in steps of delta
             double t;
             for (t = t1; t < t2 ; t += delta) X = takeStep(X, t, delta);
             if (t2 - t > 0.0) X = takeStep(X, t, t2 - t);
             return X;
         }
    }

    /** The differential equation corresponding to a trajectory with drag proportional to the square of the speed.
     * We use the fields of the Vec4 as (w, x) for position and (y, z) for velocity. */
    public static class Trajectory implements Fn {
        private final static double g = 10.0;
        private final double mass, surfaceArea, dragCoeff;

        /** Generate a differential equation giving the trajectory of an object with
         * the specified mass (from which we calculate a surface area for purposes
         * of evaluating drag) and drag coefficient. */
        public Trajectory(double _mass, double _dragCoeff) {
            mass = _mass;
            surfaceArea = Math.pow(_mass, 2.0 / 3);
            dragCoeff = _dragCoeff;
        }

        /** We take the drag force to be proportional to the surface area (and hence
         * to the mass raised to 2/3) times the velocity times the drag coefficient.
         * This gives us: dw/dt = y; dx/dt = z; dy/dt = -surfaceArea * y * dragCoeff / mass;
         * dz/dt = -surfaceArea * z * dragCoeff / mass - g */
        public Vec4 eval(double t, Vec4 X) {
           double w = X.getY();
           double x = X.getZ();
           double y = -surfaceArea * X.getY() * dragCoeff / mass;
           double z = -surfaceArea * X.getZ() * dragCoeff / mass - g;
           return new Vec4(w, x, y, z);
        }
    }

    /** Test RK4 on Trajectory by checking a known trajectory. */
    public static void main(String[] args) {
       Trajectory traj = new Trajectory(1, 0.0);
       Vec4 initialX = new Vec4(0.0, 0.0, 1.0, 10.0);
       // Should take two seconds to hit ground again
       Vec4.RK4 rk4 = new Vec4.RK4(traj);
       Vec4 finalX = rk4.integrate(initialX, 0.0, 2.0, 0.01);
       System.out.println("finalX = " + finalX);
    }
}
#9

Thank you for the source, Pjt33! I’ll try to use it tomorrow…
Concerning Semi-implicit Euler: What is the difference between it and normal Euler? I only had numerical integration mentioned to me in high school, so I can’t really understand the mathy stuff on Wiki… T___T

#10

Update the velocity using the gravitational forces at the current position. Then update the position using the new value of the velocity.

I’d be curious to know whether that makes any noticeable difference to your orbits compared to Euler with the same time step.

Doing a bit more reading, Verlet integration has the same energy-conservation property as semi-implicit Euler, and is of higher order, so that’s another thing to try.

Simon

#11

So semi-implicit Euler is just a fancy name for this code?

double dx = SUN_X - x, dy = SUN_Y - y;
double distSqrd = dx*dx + dy*dy;
double a = GRAVITATIONAL_CONSTANT * SUN_MASS / (distSqrd);
double dist = Math.sqrt(distSqrd);

velX += dx * a / dist * TIMESTEP;
velY += dy * a / dist * TIMESTEP;

x += velX * TIMESTEP;
y += velY * TIMESTEP;

That’s what I’m doing at the moment… Is normal Euler updating the position before the velocity?

#12

Yes, that looks right. If you switch the “velX/Y += …” lines with the “x/y += …” lines then you’d have standard Euler.
Does it make any difference?
Simon

#13

Made some test code to allow me to easily test things faster. So far I’ve implemented Euler, Semi-implicit Euler and Verlet (the Velocity Verlet from Wikipedia).

  • Euler: Sucks ballz. Explodes with somewhat high acceleration. Worthless for anything but exactly round orbits.
  • Semi-implicit Euler: Much better than standard Euler. Stable orbits, but the whole orbit is rotating. (see my last post)
  • Verlet: The best one by a small margin. Marginally more precise orbit compared to Semi-implicit Euler, but the exact same rotation problem. Handles extremes slightly better (only explodes at higher accelerations).

Verlet and Semi-implicit Euler followed each other almost pixel exactly for over 100k iterations, so I assume there is no significant error buildup, except for extremes where Verlet marginally survives and Euler explodes. Not really that relevant as that would just be a collision anyways, but as their computation cost is about the same, there is no reason to not use Verlet.

I still have the problem with rotating orbits. This is the most important problem to solve at the moment. I will try to implement RK4. I’ll be back, either for help or to show my results.

#14

For some reason I can’t seem to wrap my mind around RK4 for my problem. As I don’t have a function depending on time, I don’t see how I should go around implementing it… T___T I realize I should just integrate x and y separately, but I can’t seem to connect f(t) to anything decent for my problem. Gah, so confused.

#15

I did consider refactoring to remove the time parameter from the code, but I decided against. You can just ignore it.

dx/dt = vx
dy/dt = vy
dvx/dt = ax
dvy/dt = ay

You must be using those for the other quadrature implementations, so it’s just a case of putting them in a Vec4 and returning it from the eval method.

Edit: in fact, you’ve posted code further up. So

class MyFn implements Fn {
    public Vec4 eval(double t, Vec4 pos) {
        double dx = SUN_X - pos.w, dy = SUN_Y - pos.x;
        double distSqrd = dx*dx + dy*dy;
        double a = GRAVITATIONAL_CONSTANT * SUN_MASS / (distSqrd);
        double dist = Math.sqrt(distSqrd);
        return new Vec4(pos.y, pos.z, dx * a / dist, dy * a / dist);
    }
}
#16

Currently I have a Body class that has a position, velocity and a stored acceleration (needed for Verlet at least). I pass this body to an integrator using the method integrate(Body b, double timestep), which using the function getAcceleration(double x, double y) updates the body’s data. Your eval would be the same to my getAcceleration, and integrate would be takeStep, correct?
What I don’t understand in your code is how you use your Vec4s. For example:

double dx = SUN_X - pos.w, dy = SUN_Y - pos.x;

What’s up with pos.w and pos.x?

#17

I thought it would be educational (for me if no one else) to look at how some different integration methods perform for different step sizes.

So I set up a simple simulation of an asteroid following an elliptical orbit around a star, and ran it for a range of different sizes of time step. The simulations cover the length of time that the asteroid would take to complete one period of the orbit in the exact solution. So, if a simulation were perfectly accurate, the asteroid would end up exactly back at its starting point. I’m taking the distance of the asteroid from its starting point at the end of the simulation as the measure of a simulation’s error.

The asteroid starts at distance 1.0 from the star, and that’s also its maximum distance (for the exact orbit). So an error of 0.01 at the end of the simulation would probably be acceptable, although of course we’d like to do better than that. If the error is greater than 1 then the simulation was a complete mess.

The simulations run for different numbers of time steps. If T is the time simulated (i.e., one period of the exact orbit) and the simulation performs N time steps, then the size of each time step is dt=T/N.

Here are some plots showing the final error when different numbers of time steps are used. The plots cover the Euler method, the semi-implicit Euler method, and the fourth-order Runge-Kutta method. I also tried the basic Verlet and velocity Verlet methods, but the results were so similar to those of the semi-implicit Euler method that they’re not worth showing.

http://www.dishmoth.com/fig1.png

http://www.dishmoth.com/fig2.png

http://www.dishmoth.com/fig3.png

Clearly the Euler method needs a very large number of very small time steps in order for the error to be small. The semi-implicit Euler method gets away with a lot fewer time steps, and the Runge-Kutta method fewer still.

The difference in performance is clear when the error graphs are plotted all together with logarithmic axes.

http://www.dishmoth.com/fig4.png

For this plot it’s interesting to note that, when the number of time steps gets large enough, the slope of the Runge-Kutta curve is roughly -4, the slope of the semi-implicit Euler curve is roughly -2, and the slope of the Euler curve is roughly -1. This is evidence that the methods converge at fourth order, second order and first order respectively. Curiously, the semi-implicit Euler method is, like the Euler method, only supposed to converge at first order (whereas the Verlet methods converge at second order) so for this problem at least it’s performing better than advertised.

In conclusion, the Euler method isn’t worth bothering with for this problem. The semi-implicit Euler method (which is almost identical) is much better, but in terms of accuracy the Runge-Kutta method wins. That said, it takes a few times more work to perform each step of the Runge-Kutta method, and it’s harder to code, so unless you’re chasing very high accuracy it may not be worth the effort. (In effect, the green curve in the logarithmic plot should be shifted a bit to the right to account for the extra computational work involved.)

All this is just for one example orbit of course. The number of time steps needed to get decent accuracy will be different for different orbits, and will probably depend quite a lot on how close the orbit gets to the star.

I do some strange things for fun sometimes…
Simon

#18

The Vec4’s components are called w,x,y,z for generality. I’m using the same convention as in the Trajectory example class that position is (w,x) and velocity is (y,z). If you want to refactor w,x,y,z to x,y,vx,vy, feel free.

eval isn’t quite the same as getAcceleration. It gets the derivatives of all the components, which when the components are speed and velocity would be velocity and acceleration.

Nice work.

As a point of interest, my FunOrb game Torquing uses RK4 on a Vec13. IIRC the 13 components are 3 for position, 3 for velocity, 4 for angular position (using quaternions, and renormalising them regularly), and 3 for angular velocity (and I can’t remember why that isn’t using quaternions).

#19

It’s not a proper plug if you don’t give a link. :wink:

#20

It was a plug for RK4 rather than Torquing, but if you insist… :wink:
http://www.funorb.com/info.ws?game=torquing