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...

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