Help with an orbit simulator

For the discussion of math. Duh.

Moderators: gmalivuk, Moderators General, Prelates

User avatar
howardh
Posts: 100
Joined: Tue Aug 26, 2008 4:59 pm UTC

Help with an orbit simulator

Postby howardh » Sun Dec 21, 2008 4:24 pm UTC

I'm writing an orbit simulator and the goal of this project is to use the simulator to calculate the shortest path from earth to 99942 Apophis. Right now, I'm working on getting the orbits right. I'm sure they're fine when simulating each planet on their own but if I simulate two simultaneously, I think they each start at different times and so one ends up several years ahead of the other. I think this is what the epoch date is for but I'm not sure how to use it.


Also, can someone check if I applied the angles correctly?

(r, v) are the polar coordinates of the orbiting object calculated

x = cos(v + longitude of ascending node)*cos(inclination)*r
y = sin(v + longitude of ascending node)*cos(argument of perihelion)*r
z = sin(inclination)*sin(argument of perihelion)*r


When calculating the mean anomaly with the formula
(2*PI*t)/period
would using the value from this formula give the position of the object t time after the epoch with t being the same unit of time as period?


All of my information was taken off wikipedia:
http://en.wikipedia.org/wiki/99942_Apophis
http://en.wikipedia.org/wiki/Earth
http://en.wikipedia.org/wiki/Kepler%27s ... ary_motion

Here's the relevant code if anyone wants to see it:

planet.cpp

Code: Select all

void C_Planet::updatePosition(double time)
{
    t = time-this->epoch;

    l = semiLatusRectum();  //calculates the semi-latus rectum
    m = meanAnomaly();      //calculates the mean anomaly
    E = eccentricAnomaly(); //calculates the eccentric anomaly (accurate approximation using iteration)
    v = trueAnomaly();      //calculates true anomaly (for use in polar coordinates)
    r = heliocentricDistance();//calculates heliocentric distance (distance from Celestial Body, to its Gravitational centre
                                     //(body which it is in orbit around), usually the sun, hence the term 'helio'-centric distance

    r *= 20.0;

    positionPC.v = v;
    positionPC.r = r;

    positionV.x = (cos(v+loan))*c_incl*r;
    positionV.y = (sin(v+loan))*c_aop*r;
    positionV.z = s_incl*s_aop*r;
}

double C_Planet::eccentricAnomaly()
{
    // Based on code from Swift by Hal Levison and Martin Duncan
    // but they probably got it from someone else, maybe Danby.
    //assert(m>=0.);
    //assert(m<=2.*PI);
    double k = 0.85;
    double x = (m<PI) ? m + k*e : m - k*e;
    double del = 1.e-5;
    for(int i=0;i<10;++i)
    {
        double es = e*sin(x);
        double F = x-es-m;
        if(fabs(F)<del) break;
        double ec = e*cos(x);
        double Fp = 1.-ec;
        double Fpp = es;
        double Fppp = ec;
        double Dx = -F/Fp;
        Dx = -F/(Fp+Dx*Fpp/2.);
        Dx = -F/(Fp+Dx*Fpp/2.+Dx*Dx*Fppp/6.);
        x += Dx;
    }
    return x;
}
/*
double C_Planet::eccentricAnomaly()
{
    double delta = pow(10, -dp); //number of decimal places
    if (e < 0.8) {E = m;}//determines which initial value to use for E
    else         {E = PI;}//uses the radian value for 180 degrees if the eccentricity of the orbit is very high (above 0.8, which is rarely seen)
    F = E - e*sin(m) - m;//calculates the initial value to be used for iteration (function)
    while ((abs(F) > delta) && (i < maxIter)) //this for loop performs the iterations (using Newton's Method) on the Eccentric Anomaly, 'E', in Kepler's Equation
    {
            E = E - F/(1.0 - e*cos(E));//iteration of eccentric anomaly
            F = E - e*sin(E) - m;  //calculates the value to be used for the next iteration
            i++;    //incriments counter
    }
    //return (E * pow(10, dp)) / pow(10, dp); //rounds to the requested number of decimal places
    return E;
}
*/
double C_Planet::meanAnomaly() //calculates the mean anomaly (see Kepler's Equation)
{
    return (2*(PI)*t) / p;
}

double C_Planet::trueAnomaly() //calculates true anomly (see Kepler's Equation) for use in polar coordinates
{
    return (atan(sqrt((1 + e)/(1 - e))*tan(E/2)))*2;
}

double C_Planet::heliocentricDistance() //calculates heliocentric distance (see Kepler's Equation) for use in polar coordinates
{
    return a*(1-e*e)/(1+e*cos(v));
}

double C_Planet::semiLatusRectum() //calculates semi-latus rectum (distance from centre of ellipse to the focii), used to calculate the heliocentric distance
{
    return a*(1 - (e*e));
}

void C_Planet::setInclination(double d)
{
    this->incl = d;
    this->s_incl = sin(d);
    this->c_incl = cos(d);
}

void C_Planet::setLoan(double d)
{
    this->loan = d;
    this->s_loan = sin(d);
    this->c_loan = cos(d);
}

void C_Planet::setAop(double d)
{
    this->aop = d;
    this->s_aop = sin(d);
    this->c_aop = cos(d);
}


planet.h

Code: Select all

class C_Planet
{
    public:
        C_Planet():dp(5),i(0),maxIter(30) {}
        C_Planet(double ecc, double per, double sm):dp(5),i(0),maxIter(30), e(ecc), p(per), a(sm) {}

        void updatePosition(double time);
        CVector3 getPosition() { return positionV; }
        PolarCoord getPositionPC() { return positionPC; }

        double eccentricAnomaly();
        double meanAnomaly();
        double trueAnomaly();
        double heliocentricDistance();
        double semiLatusRectum();

        void setEpoch(double d) { this->epoch = d; }

        void setInclination(double d);
        void setLoan(double d);
        void setAop(double d);

        double getSinIncl() { return s_incl; }
        double getSinLoan() { return s_loan; }
        double getSinAop()  { return s_aop; }

        double getCosIncl() { return c_incl; }
        double getCosLoan() { return c_loan; }
        double getCosAop()  { return c_aop; }

        void setRadius(double rad) { radius = rad; }
        double getRadius() { return radius; }
    private:
        CVector3 positionV;    //Cartesian coordinates of the planet position
        PolarCoord positionPC; //Polar coordinates of the planet position

        double epoch; //epoch date in Julian Days

        double m;    //mean anomaly (for all terms see wikipedia page on 'Kepler's Laws'
        double e;    //eccentricity of orbit
        double E;    //Eccentric anomaly
        double F;    //(function), used to store value during iteration
        double p;    //period length
        double t;    //time since perihelion
        double r;    //heliocentric distance (distance from Sun to Body, in Astronomical Units, used for polar coordinate grid)
        double v;    //angle (true anomaly) centred at the Sun, from the point of perihelion (closest pt. to Sun) to location of Body, calculated in radians
        double a;    //semi-major axis of ellipse
        double l;    //semi-latus rectum (dist from centre to foci)
        int dp;      //# of decimal places for iteration
        int i;       //counter for iterations
        int maxIter; //max number of iterations

        double incl; //inclination: rotate on y acis
        double loan; //longitude of ascending node: rotate on z axis
        double aop;  //argument of perihelion: rotate on x axis

        //sine
        double s_incl; //inclination
        double s_loan; //longitude of ascending node
        double s_aop;  //argument of perihelion

        //cosine
        double c_incl; //inclination
        double c_loan; //longitude of ascending node
        double c_aop;  //argument of perihelion

        double radius; //radius in AU * 1000
};


edit:
I'm saddened by the lack of replies... :cry:
A small change was made to the code (epoch) though I'm not entirely sure that that part's done right either
I tested the simulation several times but never got the expected results. I checked the code converting the polar coordinates to cartesian corrdinates by comparing the magnitude after the conversion and the r value of the polar coordinate, the value changed so I did mess up there.

edit2:

Code: Select all

    positionV.x = (cos(v+loan))*c_incl*r;
    positionV.y = (sin(v+loan))*r;
    positionV.z = -(cos(v+loan))*s_incl*r;

    positionV.y = positionV.y*c_aop-positionV.z*s_aop;
    positionV.z = positionV.z*c_aop+positionV.y*s_aop;

Still not getting the expected results
meow

User avatar
scikidus
Posts: 792
Joined: Wed Nov 26, 2008 9:34 pm UTC
Location: New York, NY
Contact:

Re: Help with an orbit simulator

Postby scikidus » Thu Dec 25, 2008 3:11 pm UTC

I'm not entirely sure if I can help at all, but a few weeks ago I built a Newtonian Physics simulator similar to this. All it uses is the formula F=GMm/r^2 and calculates the respective forces between all objects to move them accordingly.

http://freewebs.com/scikidus/NewPhy.swf

You can throw objects (hit space to start up the simulator again), and resize objects by holding down Control and clicking.

If you're interested, I can post up the code. (It's in Flash.)
Happy hollandaise!

"The universe is a figment of its own imagination" -Douglas Adams

User avatar
howardh
Posts: 100
Joined: Tue Aug 26, 2008 4:59 pm UTC

Re: Help with an orbit simulator

Postby howardh » Fri Dec 26, 2008 1:04 pm UTC

Ooh, a reply...
That wouldn't help much here but thanks anyways.
meow

jr_naxera
Posts: 2
Joined: Wed Jan 14, 2009 7:47 pm UTC

Re: Help with an orbit simulator

Postby jr_naxera » Wed Jan 14, 2009 7:58 pm UTC

Hey howardh, I'm working on a similar project in C++, only mine is oriented to eventually became the physics engine for a strategy game, so precision is not much of a concern for me. However, I have found some interesting links, this one in particular seems suited to serve your needs.

http://jenab6.livejournal.com/12053.html

There is also a page with excellent formulas for orbital calculations, I'm gonna search it later and post it.

Cheers.

jr_naxera
Posts: 2
Joined: Wed Jan 14, 2009 7:47 pm UTC

Re: Help with an orbit simulator

Postby jr_naxera » Thu Jan 15, 2009 2:55 am UTC

Well, as promised, this page has every formula you may need to describe the orbit from the six classical keplerian elements. Just be careful because it uses a slightly confusing nomenclature:

http://www.nature1st.net/bogan/orbits/kepler/orbteqtn.html

Also these are my routines to translate from the planar polar coordinates of the orbit to heliocentric coordinates, I took them from a very useful text in Google Books, called "Fundamentals of Astrodynamics and Aplications" by David Vallado, where the math is really well explained:

Code: Select all

//This function calculate the conversion factors using euler angles
//alpha is the Right Ascencion of the Ascending Node (Raan)
//beta is the Inclination
//gamma is the Argument of Periapsis
//euler_f is a structure for passing a rotation matrix
euler_f rot_factors(double alpha, double beta, double gamma)
{
    euler_f t;
    t.x1=  (cos(alpha)*cos(gamma))
          -(sin(alpha)*sin(gamma)*cos(beta));
    t.x2= -(cos(alpha)*sin(gamma))
          -(sin(alpha)*cos(gamma)*cos(beta));
    t.y1=  (sin(alpha)*cos(gamma))
          +(cos(alpha)*sin(gamma)*cos(beta));
    t.y2= -(sin(alpha)*sin(gamma))
          +(cos(alpha)*cos(gamma)*cos(beta));
    t.z1= sin(gamma)*sin(beta);
    t.z2= cos(gamma)*sin(beta);

    return t;
}
//This function calculate the position vector using the previous factors and the polar coordinates from the keplerian elements
//k is a structure for passing the keplerian elements, with k.dperi as the periapsis distance and k.anomaly as the true anomaly
//eu_vector is an object that represents an Euclidean Vector (used to abstract all vector-related operations and conversions)
eu_vector kepler_to_r(keplerians &k, euler_f &f)
{
    double fr1=(k.dperi*cos(k.anomaly))/(1+(k.e*cos(k.anomaly)));
    double fr2=(k.dperi*sin(k.anomaly))/(1+(k.e*cos(k.anomaly)));
    double rx=(fr1*f.x1)+(fr2*f.x2);
    double ry=(fr1*f.y1)+(fr2*f.y2);
    double rz=(fr1*f.z1)+(fr2*f.z2);

    return eu_vector(rx,ry,rz,CART);
}
//This function calculate the velocity vector using the factors and the polar coordinates from the keplerian elements
//focus is an object that contains the information of the point mass around which the orbit is traced, with Mu() as its Standard Gravitational Parameter
eu_vector kepler_to_v(keplerians &k, point_mass &focus, euler_f &f)//Función para obtener el Vector de Velocidad
{
    double fv1=-(sqrt(focus.Mu()/k.dperi)*sin(k.anomaly));
    double fv2=  sqrt(focus.Mu()/k.dperi)*(k.e+cos(k.anomaly));
    double vx=(fv1*f.x1)+(fv2*f.x2);
    double vy=(fv1*f.y1)+(fv2*f.y2);
    double vz=(fv1*f.z1)+(fv2*f.z2);

    return eu_vector(vx,vy,vz,CART);
}


I hope you find some of this information useful, it surely has been a great help to me.

Cheers

User avatar
howardh
Posts: 100
Joined: Tue Aug 26, 2008 4:59 pm UTC

Re: Help with an orbit simulator

Postby howardh » Mon Jan 19, 2009 11:06 pm UTC

I've done away with Kepler's equation. The two bodies kept going into opposite ends of the Solar System when they should be colliding, probably because of inaccuracies in using this method, or maybe I just didn't implement it properly, I don't know and it doesn't matter anymore.
I have the simulation working properly using NASA's SPICE toolkit and plan on using F=G*(m1*m2)/r^2 like in scikidus' simulator for the next part of the project.

Thanks for your replies.

I'll be back once I hit another road block...
meow

User avatar
ConMan
Shepherd's Pie?
Posts: 1690
Joined: Tue Jan 01, 2008 11:56 am UTC
Location: Beacon Alpha

Re: Help with an orbit simulator

Postby ConMan » Tue Jan 20, 2009 4:41 am UTC

howardh wrote:I've done away with Kepler's equation. The two bodies kept going into opposite ends of the Solar System when they should be colliding, probably because of inaccuracies in using this method, or maybe I just didn't implement it properly, I don't know and it doesn't matter anymore.
I have the simulation working properly using NASA's SPICE toolkit and plan on using F=G*(m1*m2)/r^2 like in scikidus' simulator for the next part of the project.

Thanks for your replies.

I'll be back once I hit another road block...


It might be a rounding error effect. Euler's Method performs particularly poorly on things like orbits because it doesn't conserve what are typically considered to be conserved quantities (like angular momentum). Try reading up on the Verlet and Runje-Kutta methods to see if they prove more stable.
pollywog wrote:
Wikihow wrote:* Smile a lot! Give a gay girl a knowing "Hey, I'm a lesbian too!" smile.
I want to learn this smile, perfect it, and then go around smiling at lesbians and freaking them out.

stephentyrone
Posts: 778
Joined: Mon Aug 11, 2008 10:58 pm UTC
Location: Palo Alto, CA

Re: Help with an orbit simulator

Postby stephentyrone » Tue Jan 20, 2009 5:31 am UTC

ConMan wrote:Try reading up on the Verlet and Runje-Kutta methods to see if they prove more stable.


The key phrase for talking about solvers that conserve certain quantities is "symplectic methods". To find out more, see any of a number of books on Geometric Numerical Integration, several papers by William Kahan and others, and try the midpoint rule.

Also some runge-kutta methods are symplectic, but most are not.
GENERATION -16 + 31i: The first time you see this, copy it into your sig on any forum. Square it, and then add i to the generation.

User avatar
taby
Posts: 154
Joined: Thu Nov 01, 2007 8:39 pm UTC
Location: Canada

Re: Help with an orbit simulator

Postby taby » Tue Jan 20, 2009 2:47 pm UTC

re: Runge-Kutta
http://gafferongames.wordpress.com/game ... on-basics/
http://gafferongames.wordpress.com/game ... -timestep/

If you get stuck, I can try to help. I've got source for this type of example handy in C++ (CPU) and GLSL (GPU) using Euler, RK2, and RK4 integration, if you need it. It includes a simple camera class, all in OpenGL of course.

I'm not sure if this is obvious or not, but one will find that the behaviour described by Kepler's equations matches the behaviour described by Newton's F = GMm/r^2.


Return to “Mathematics”

Who is online

Users browsing this forum: No registered users and 6 guests