Orbit.cpp   Orbit.cpp 
skipping to change at line 45 skipping to change at line 45
#define cbrt(x) pow((x),1./3.) #define cbrt(x) pow((x),1./3.)
#endif #endif
//! Solve true anomaly nu for hyperbolic orbit. //! Solve true anomaly nu for hyperbolic orbit.
//! @param q: perihel distance //! @param q: perihel distance
//! @param n: mean motion //! @param n: mean motion
//! @param e: excentricity //! @param e: excentricity
//! @param dt: days from perihel //! @param dt: days from perihel
//! @param rCosNu: r*cos(nu) //! @param rCosNu: r*cos(nu)
//! @param rSinNu: r*sin(nu) //! @param rSinNu: r*sin(nu)
static void InitHyp(const double q, const double n, const double e, const d static void InitHyp(const double q, const double n, const double e, const d
ouble dt, double &rCosNu, double &rSinNu) { ouble dt, double &rCosNu, double &rSinNu)
// qDebug() << "InitHyp"; {
Q_ASSERT(e>1.0); // qDebug() << "InitHyp";
const double a = q/(e-1.0); Q_ASSERT(e>1.0);
Q_ASSERT(a>0.0); const double a = q/(e-1.0);
const double M = n * dt; Q_ASSERT(a>0.0);
// double H = M; const double M = n * dt;
// for (;;) { // Newton // double H = M;
// const double Hp = H; // for (;;) { // Newton
// H = H-(e*sinh(H)-H-M)/(e*cosh(H)-1); // const double Hp = H;
// if (fabs(H - Hp) < EPSILON) break; // H = H-(e*sinh(H)-H-M)/(e*cosh(H)-1);
// } // if (fabs(H - Hp) < EPSILON) break;
// const double h1 = q*sqrt((e+1.0)/(e-1.0)); // }
// a1 = a*(e-cosh(H)); // const double h1 = q*sqrt((e+1.0)/(e-1.0));
// a2 = h1*sinh(H); // a1 = a*(e-cosh(H));
// GZ Again I prefer Heafner, ch.5.4 // a2 = h1*sinh(H);
double E=StelUtils::sign(M)*log(2.0*fabs(M)/e + 1.85); // GZ Again I prefer Heafner, ch.5.4
// qDebug() << "InitHyp: E=" << E << " M=" << M ; double E=StelUtils::sign(M)*log(2.0*fabs(M)/e + 1.85);
for (;;){ // qDebug() << "InitHyp: E=" << E << " M=" << M ;
const double Ep=E; for (;;)
const double f2=e*sinh(E); {
const double f=f2-E-M; const double Ep=E;
const double f1=e*cosh(E)-1.0; const double f2=e*sinh(E);
E+= (-5.0*f)/(f1+StelUtils::sign(f1)*sqrt(fabs(16.0*f1*f1-20.0*f*f2)) const double f=f2-E-M;
); const double f1=e*cosh(E)-1.0;
if (fabs(E-Ep) < EPSILON) break; E+= (-5.0*f)/(f1+StelUtils::sign(f1)*std::sqrt(fabs(16.0*f1*
} f1-20.0*f*f2)));
rCosNu = a*(e-cosh(E)); if (fabs(E-Ep) < EPSILON) break;
rSinNu = a*sqrt(e*e-1.0)*sinh(E); }
rCosNu = a*(e-cosh(E));
rSinNu = a*std::sqrt(e*e-1.0)*sinh(E);
} }
//! Solve true anomaly nu for parabolic orbit. //! Solve true anomaly nu for parabolic orbit.
//! @param q: perihel distance //! @param q: perihel distance
//! @param n: mean motion equivalent related to W (n=W/dt) in Heafner, ch5. 5 //! @param n: mean motion equivalent related to W (n=W/dt) in Heafner, ch5. 5
//! @param dt: days from perihel //! @param dt: days from perihel
//! @param rCosNu: r*cos(nu) //! @param rCosNu: r*cos(nu)
//! @param rSinNu: r*sin(nu) //! @param rSinNu: r*sin(nu)
/* /*
static void InitPar(const double q, const double n, const double dt, doub static void InitPar(const double q, const double n, const double dt, double
le &a1, double &a2) { &a1, double &a2)
const double A = n*dt; {
const double h = sqrt(A*A+1.0); const double A = n*dt;
double c = cbrt(fabs(A)+h); const double h = sqrt(A*A+1.0);
c = c*c; double c = cbrt(fabs(A)+h);
const double tan_nu_h = 2*A/(1+c+1/c); c = c*c;
a1 = q*(1-tan_nu_h*tan_nu_h); const double tan_nu_h = 2*A/(1+c+1/c);
a2 = 2.0*q*tan_nu_h; a1 = q*(1-tan_nu_h*tan_nu_h);
a2 = 2.0*q*tan_nu_h;
} }
*/ */
// GZ This implementation now follows Heafner, ch 5.5 // GZ This implementation now follows Heafner, ch 5.5
static void InitPar(const double q, const double n, const double dt, double static void InitPar(const double q, const double n, const double dt, double
&rCosNu, double &rSinNu) { &rCosNu, double &rSinNu)
// qDebug() << "InitPar"; {
// const double M=dt*sqrt(GAUSS_GRAV_CONST/(2.0*q*q*q)); // qDebug() << "InitPar";
// const double W=1.5*M; // const double M=dt*sqrt(GAUSS_GRAV_CONST/(2.0*q*q*q));
const double W=dt*n; // const double W=1.5*M;
const double Y=cbrt(W+sqrt(W*W+1)); const double W=dt*n;
const double tanNu2=Y-1.0/Y; // Heafner (5.5.8) has an error here, writ const double Y=cbrt(W+std::sqrt(W*W+1));
es (Y-1)/Y. const double tanNu2=Y-1.0/Y; // Heafner (5.5.8) has an error here, w
rCosNu=q*(1.0-tanNu2*tanNu2); rites (Y-1)/Y.
rSinNu=2.0*q*tanNu2; rCosNu=q*(1.0-tanNu2*tanNu2);
rSinNu=2.0*q*tanNu2;
} }
//! Solve true anomaly nu for elliptical orbit with Laguerre-Conway's metho d. (May have high e) //! Solve true anomaly nu for elliptical orbit with Laguerre-Conway's metho d. (May have high e)
//! @param q: perihel distance //! @param q: perihel distance
//! @param n: mean motion //! @param n: mean motion
//! @param e: excentricity //! @param e: excentricity
//! @param dt: days from perihel //! @param dt: days from perihel
//! @param rCosNu: r*cos(nu) //! @param rCosNu: r*cos(nu)
//! @param rSinNu: r*sin(nu) //! @param rSinNu: r*sin(nu)
static void InitEll(const double q, const double n, const double e, const d static void InitEll(const double q, const double n, const double e, const d
ouble dt, double &rCosNu, double &rSinNu) { ouble dt, double &rCosNu, double &rSinNu)
// qDebug() << "InitEll"; {
Q_ASSERT(e<1.0); // qDebug() << "InitEll";
const double a = q/(1.0-e); // semimajor axis Q_ASSERT(e<1.0);
double M = fmod(n*dt,2*M_PI); // Mean Anomaly const double a = q/(1.0-e); // semimajor axis
if (M < 0.0) M += 2.0*M_PI; double M = fmod(n*dt,2*M_PI); // Mean Anomaly
// double E = M; if (M < 0.0) M += 2.0*M_PI;
// for (;;) { // Newton(?) Solve Kepler's equation (similar to Meeus secon // double E = M;
d method, Astro.Alg. 1998 p.199) // for (;;) { // Newton(?) Solve Kepler's equation (similar to Meeus se
// const double Ep = E; cond method, Astro.Alg. 1998 p.199)
// E -= (M-E+e*sin(E))/(e*cos(E)-1); // const double Ep = E;
// if (fabs(E-Ep) < EPSILON) break; // E -= (M-E+e*sin(E))/(e*cos(E)-1);
// } // if (fabs(E-Ep) < EPSILON) break;
// GZ: Comet orbits are quite often near-parabolic, where this may still // }
only converge slowly. // GZ: Comet orbits are quite often near-parabolic, where this may stil
// Better always use Laguerre-Conway. See Heafner, Ch. 5.3 l only converge slowly.
double E=M+0.85*e*StelUtils::sign(sin(M)); // Better always use Laguerre-Conway. See Heafner, Ch. 5.3
for (;;){ double E=M+0.85*e*StelUtils::sign(sin(M));
const double Ep=E; for (;;)
const double f2=e*sin(E); {
const double f=E-f2-M; const double Ep=E;
const double f1=1.0-e*cos(E); const double f2=e*sin(E);
E+= (-5.0*f)/(f1+StelUtils::sign(f1)*sqrt(fabs(16.0*f1*f1-20.0*f*f2)) const double f=E-f2-M;
); const double f1=1.0-e*cos(E);
if (fabs(E-Ep) < EPSILON) break; E+= (-5.0*f)/(f1+StelUtils::sign(f1)*std::sqrt(fabs(16.0*f1*
} f1-20.0*f*f2)));
// Note: q=a*(1-e) if (fabs(E-Ep) < EPSILON) break;
const double h1 = q*sqrt((1.0+e)/(1.0-e)); // elsewhere: a sqrt(1-e²) }
... q / (1-e) sqrt( (1+e)(1-e)) = q sqrt((1+e)/(1-e)) // Note: q=a*(1-e)
rCosNu = a*(cos(E)-e); const double h1 = q*std::sqrt((1.0+e)/(1.0-e)); // elsewhere: a sqr
rSinNu = h1*sin(E); t(1-e²) ... q / (1-e) sqrt( (1+e)(1-e)) = q sqrt((1+e)/(1-e))
rCosNu = a*(cos(E)-e);
rSinNu = h1*sin(E);
} }
//! Compute position vector and (optional) speed vector from orbital elemen ts and true anomaly components. See e.g. Heafner, Fund.Eph.Comp.1999 //! Compute position vector and (optional) speed vector from orbital elemen ts and true anomaly components. See e.g. Heafner, Fund.Eph.Comp.1999
//! @param i inclination //! @param i inclination
//! @param Omega, longitude of ascending node //! @param Omega, longitude of ascending node
//! @param w omega, argument of pericenter //! @param w omega, argument of pericenter
//! @param rCosNu: r*cos(nu) //! @param rCosNu: r*cos(nu)
//! @param rSinNu: r*sin(nu) //! @param rSinNu: r*sin(nu)
//! @param rx: x component of position vector, AU //! @param rx: x component of position vector, AU
//! @param ry: y component of position vector, AU //! @param ry: y component of position vector, AU
//! @param rz: z component of position vector, AU //! @param rz: z component of position vector, AU
//! @param withVelVector also compute velocity vector (required for comet t ails) //! @param withVelVector also compute velocity vector (required for comet t ails)
//! @param e excentricity (required if withVelVector=true) //! @param e excentricity (required if withVelVector=true)
//! @param q perihel distance, AU (required if withVelVector=true) //! @param q perihel distance, AU (required if withVelVector=true)
//! @param rdotx: x component of velocity vector, AU/d //! @param rdotx: x component of velocity vector, AU/d
//! @param rdoty: y component of velocity vector, AU/d //! @param rdoty: y component of velocity vector, AU/d
//! @param rdotz: z component of velocity vector, AU/d //! @param rdotz: z component of velocity vector, AU/d
void Init3D(const double i, const double Omega, const double w, const doubl e rCosNu, const double rSinNu, void Init3D(const double i, const double Omega, const double w, const doubl e rCosNu, const double rSinNu,
double &rx,double &ry,double &rz, double &rdotx, double &rdoty, double &rx,double &ry,double &rz, double &rdotx, double &rdoty,
double &rdotz, const bool withVelVector=false, const double e=0.0, const d double &rdotz, const bool withVelVector=false,
ouble q=0.0) { const double e=0.0, const double q=0.0)
const double cw = cos(w); {
const double sw = sin(w); const double cw = cos(w);
const double cOm = cos(Omega); const double sw = sin(w);
const double sOm = sin(Omega); const double cOm = cos(Omega);
const double ci = cos(i); const double sOm = sin(Omega);
const double si = sin(i); const double ci = cos(i);
const double Px=-sw*sOm*ci+cw*cOm; // Heafner, 5.3.1 Px const double si = sin(i);
const double Qx=-cw*sOm*ci-sw*cOm; // Heafner, 5.3.4 Qx const double Px=-sw*sOm*ci+cw*cOm; // Heafner, 5.3.1 Px
const double Py= sw*cOm*ci+cw*sOm; // Heafner, 5.3.2 Py const double Qx=-cw*sOm*ci-sw*cOm; // Heafner, 5.3.4 Qx
const double Qy= cw*cOm*ci-sw*sOm; // Heafner, 5.3.5 Qy const double Py= sw*cOm*ci+cw*sOm; // Heafner, 5.3.2 Py
const double Pz= sw*si; // Heafner, 5.3.3 Pz const double Qy= cw*cOm*ci-sw*sOm; // Heafner, 5.3.5 Qy
const double Qz= cw*si; // Heafner, 5.3.6 Qz const double Pz= sw*si; // Heafner, 5.3.3 Pz
rx = Px*rCosNu+Qx*rSinNu; // Heafner, 5.3.18 r const double Qz= cw*si; // Heafner, 5.3.6 Qz
ry = Py*rCosNu+Qy*rSinNu; rx = Px*rCosNu+Qx*rSinNu; // Heafner, 5.3.18 r
rz = Pz*rCosNu+Qz*rSinNu; ry = Py*rCosNu+Qy*rSinNu;
if (withVelVector) { rz = Pz*rCosNu+Qz*rSinNu;
const double r=sqrt(rSinNu*rSinNu+rCosNu*rCosNu); if (withVelVector)
const double sinNu=rSinNu/r; {
const double cosNu=rCosNu/r; const double r=std::sqrt(rSinNu*rSinNu+rCosNu*rCosNu);
const double p=q*(1.0+e); const double sinNu=rSinNu/r;
const double sqrtMuP=sqrt(GAUSS_GRAV_CONST/p); const double cosNu=rCosNu/r;
rdotx=sqrtMuP*((e+cosNu)*Qx - sinNu*Px); // Heafner, 5.3.19 r' const double p=q*(1.0+e);
rdoty=sqrtMuP*((e+cosNu)*Qy - sinNu*Py); const double sqrtMuP=std::sqrt(GAUSS_GRAV_CONST/p);
rdotz=sqrtMuP*((e+cosNu)*Qz - sinNu*Pz); rdotx=sqrtMuP*((e+cosNu)*Qx - sinNu*Px); // Heafner, 5.3.19
} r'
rdoty=sqrtMuP*((e+cosNu)*Qy - sinNu*Py);
rdotz=sqrtMuP*((e+cosNu)*Qz - sinNu*Pz);
}
} }
CometOrbit::CometOrbit(double pericenterDistance, CometOrbit::CometOrbit(double pericenterDistance,
double eccentricity, double eccentricity,
double inclination, double inclination,
double ascendingNode, double ascendingNode,
double argOfPerhelion, double argOfPerhelion,
double timeAtPerihelion, double timeAtPerihelion,
double orbitGoodDays, double orbitGoodDays,
double meanMotion, // GZ: for parabolic s, this is W/dt in Heafner's lettering double meanMotion, // GZ: for parabolic s, this is W/dt in Heafner's lettering
double parentRotObliquity, double parentRotObliquity,
double parentRotAscendingnode, double parentRotAscendingnode,
double parentRotJ2000Longitude) double parentRotJ2000Longitude)
//) //)
:q(pericenterDistance),e(eccentricity),i(inclination), : q(pericenterDistance),
Om(ascendingNode),w(argOfPerhelion),t0(timeAtPerihelion), e(eccentricity),
n(meanMotion), updateTails(true), orbitGood(orbitGoodDays) { i(inclination),
rdot.set(0.0, 0.0, 0.0); Om(ascendingNode),
const double c_obl = cos(parentRotObliquity); w(argOfPerhelion),
const double s_obl = sin(parentRotObliquity); t0(timeAtPerihelion),
const double c_nod = cos(parentRotAscendingnode); n(meanMotion),
const double s_nod = sin(parentRotAscendingnode); updateTails(true),
const double cj = cos(parentRotJ2000Longitude); orbitGood(orbitGoodDays)
const double sj = sin(parentRotJ2000Longitude); {
rdot.set(0.0, 0.0, 0.0);
// rotateToVsop87[0] = c_nod; const double c_obl = cos(parentRotObliquity);
// rotateToVsop87[1] = -s_nod * c_obl; const double s_obl = sin(parentRotObliquity);
// rotateToVsop87[2] = s_nod * s_obl; const double c_nod = cos(parentRotAscendingnode);
// rotateToVsop87[3] = s_nod; const double s_nod = sin(parentRotAscendingnode);
// rotateToVsop87[4] = c_nod * c_obl; const double cj = cos(parentRotJ2000Longitude);
// rotateToVsop87[5] = -c_nod * s_obl; const double sj = sin(parentRotJ2000Longitude);
// rotateToVsop87[6] = 0.0;
// rotateToVsop87[7] = s_obl; // rotateToVsop87[0] = c_nod;
// rotateToVsop87[8] = c_obl; // rotateToVsop87[1] = -s_nod * c_obl;
rotateToVsop87[0] = c_nod*cj-s_nod*c_obl*sj; // rotateToVsop87[2] = s_nod * s_obl;
rotateToVsop87[1] = -c_nod*sj-s_nod*c_obl*cj; // rotateToVsop87[3] = s_nod;
rotateToVsop87[2] = s_nod*s_obl; // rotateToVsop87[4] = c_nod * c_obl;
rotateToVsop87[3] = s_nod*cj+c_nod*c_obl*sj; // rotateToVsop87[5] = -c_nod * s_obl;
rotateToVsop87[4] = -s_nod*sj+c_nod*c_obl*cj; // rotateToVsop87[6] = 0.0;
rotateToVsop87[5] = -c_nod*s_obl; // rotateToVsop87[7] = s_obl;
rotateToVsop87[6] = s_obl*sj; // rotateToVsop87[8] = c_obl;
rotateToVsop87[7] = s_obl*cj; rotateToVsop87[0] = c_nod*cj-s_nod*c_obl*sj;
rotateToVsop87[8] = c_obl; rotateToVsop87[1] = -c_nod*sj-s_nod*c_obl*cj;
// qDebug() << "CometOrbit::()...done"; rotateToVsop87[2] = s_nod*s_obl;
} rotateToVsop87[3] = s_nod*cj+c_nod*c_obl*sj;
rotateToVsop87[4] = -s_nod*sj+c_nod*c_obl*cj;
void CometOrbit::positionAtTimevInVSOP87Coordinates(double JD, double *v, b rotateToVsop87[5] = -c_nod*s_obl;
ool updateVelocityVector) { rotateToVsop87[6] = s_obl*sj;
JD -= t0; rotateToVsop87[7] = s_obl*cj;
double rCosNu,rSinNu; rotateToVsop87[8] = c_obl;
// temporary solve freezes for near-parabolic comets - using (e < 0.9999) // qDebug() << "CometOrbit::()...done";
for elliptical orbits }
// TODO: improve calculations orbits for near-parabolic comets --AW
// if (e < 0.9999) InitEll(q,n,e,JD,a1,a2); void CometOrbit::positionAtTimevInVSOP87Coordinates(double JD, double *v, b
if (e < 1.0) InitEll(q,n,e,JD,rCosNu,rSinNu); // GZ: After solving with L ool updateVelocityVector)
aguerre-Conway, I dare to go for 1.0. {
else if (e > 1.0) { JD -= t0;
// qDebug() << "Hyperbolic orbit for ecc=" << e << ", i=" << i << ", w double rCosNu,rSinNu;
=" << w << ", Mean Motion n=" << n; // temporary solve freezes for near-parabolic comets - using (e < 0.999
InitHyp(q,n,e,JD,rCosNu,rSinNu); 9) for elliptical orbits
} // TODO: improve calculations orbits for near-parabolic comets --AW
else InitPar(q,n,JD,rCosNu,rSinNu); // if (e < 0.9999) InitEll(q,n,e,JD,a1,a2);
double p0,p1,p2, s0, s1, s2; if (e < 1.0) InitEll(q,n,e,JD,rCosNu,rSinNu); // GZ: After solving w
Init3D(i,Om,w,rCosNu,rSinNu,p0,p1,p2, s0, s1, s2, updateVelocityVector, e ith Laguerre-Conway, I dare to go for 1.0.
, q); else if (e > 1.0)
v[0] = rotateToVsop87[0]*p0 + rotateToVsop87[1]*p1 + rotateToVsop87[2]*p2 {
; // qDebug() << "Hyperbolic orbit for ecc=" << e << ", i=" <<
v[1] = rotateToVsop87[3]*p0 + rotateToVsop87[4]*p1 + rotateToVsop87[5]*p2 i << ", w=" << w << ", Mean Motion n=" << n;
; InitHyp(q,n,e,JD,rCosNu,rSinNu);
v[2] = rotateToVsop87[6]*p0 + rotateToVsop87[7]*p1 + rotateToVsop87[8]*p2 }
; else InitPar(q,n,JD,rCosNu,rSinNu);
double p0,p1,p2, s0, s1, s2;
if (updateVelocityVector) { Init3D(i,Om,w,rCosNu,rSinNu,p0,p1,p2, s0, s1, s2, updateVelocityVect
rdot.set(s0, s1, s2); or, e, q);
updateTails=true; v[0] = rotateToVsop87[0]*p0 + rotateToVsop87[1]*p1 + rotateToVsop87[
} 2]*p2;
v[1] = rotateToVsop87[3]*p0 + rotateToVsop87[4]*p1 + rotateToVsop87[
5]*p2;
v[2] = rotateToVsop87[6]*p0 + rotateToVsop87[7]*p1 + rotateToVsop87[
8]*p2;
if (updateVelocityVector)
{
rdot.set(s0, s1, s2);
updateTails=true;
}
} }
EllipticalOrbit::EllipticalOrbit(double pericenterDistance, EllipticalOrbit::EllipticalOrbit(double pericenterDistance,
double eccentricity, double eccentricity,
double inclination, double inclination,
double ascendingNode, double ascendingNode,
double argOfPeriapsis, double argOfPeriapsis,
double meanAnomalyAtEpoch, double meanAnomalyAtEpoch,
double period, double period,
double epoch, double epoch,
double parentRotObliquity, double parentRotObliquity,
double parentRotAscendingnode, double parentRotAscendingnode,
double pare double parentRotJ2000Longitude)
ntRotJ2000Longitude) : : pericenterDistance(pericenterDistance),
eccentricity(eccentricity),
pericenterDistance(pericenterDistance), inclination(inclination),
eccentricity(eccentricity), ascendingNode(ascendingNode),
inclination(inclination), argOfPeriapsis(argOfPeriapsis),
ascendingNode(ascendingNode), meanAnomalyAtEpoch(meanAnomalyAtEpoch),
argOfPeriapsis(argOfPeriapsis), period(period),
meanAnomalyAtEpoch(meanAnomalyAtEpoch), epoch(epoch)
period(period), {
epoch(epoch) const double c_obl = cos(parentRotObliquity);
{ const double s_obl = sin(parentRotObliquity);
const double c_nod = cos(parentRotAscendingnode);
const double c_obl = cos(parentRotObliquity); const double s_nod = sin(parentRotAscendingnode);
const double s_obl = sin(parentRotObliquity); const double cj = cos(parentRotJ2000Longitude);
const double c_nod = cos(parentRotAscendingnode); const double sj = sin(parentRotJ2000Longitude);
const double s_nod = sin(parentRotAscendingnode);
const double cj = cos(parentRotJ2000Longitude); rotateToVsop87[0] = c_nod*cj-s_nod*c_obl*sj;
const double sj = sin(parentRotJ2000Longitude); rotateToVsop87[1] = -c_nod*sj-s_nod*c_obl*cj;
rotateToVsop87[2] = s_nod*s_obl;
rotateToVsop87[0] = c_nod*cj-s_nod*c_obl*sj; rotateToVsop87[3] = s_nod*cj+c_nod*c_obl*sj;
rotateToVsop87[1] = -c_nod*sj-s_nod*c_obl*cj; rotateToVsop87[4] = -s_nod*sj+c_nod*c_obl*cj;
rotateToVsop87[2] = s_nod*s_obl; rotateToVsop87[5] = -c_nod*s_obl;
rotateToVsop87[3] = s_nod*cj+c_nod*c_obl*sj; rotateToVsop87[6] = s_obl*sj;
rotateToVsop87[4] = -s_nod*sj+c_nod*c_obl*cj; rotateToVsop87[7] = s_obl*cj;
rotateToVsop87[5] = -c_nod*s_obl; rotateToVsop87[8] = c_obl;
rotateToVsop87[6] = s_obl*sj;
rotateToVsop87[7] = s_obl*cj;
rotateToVsop87[8] = c_obl;
} }
// Standard iteration for solving Kepler's Equation // Standard iteration for solving Kepler's Equation
struct SolveKeplerFunc1 : public unary_function<double, double> struct SolveKeplerFunc1 : public unary_function<double, double>
{ {
double ecc; double ecc;
double M; double M;
SolveKeplerFunc1(double _ecc, double _M) : ecc(_ecc), M(_M) {} SolveKeplerFunc1(double _ecc, double _M) : ecc(_ecc), M(_M) {}
double operator()(double x) const double operator()(double x) const
{ {
return M + ecc * sin(x); return M + ecc * sin(x);
} }
}; };
// Faster converging iteration for Kepler's Equation; more efficient // Faster converging iteration for Kepler's Equation; more efficient
// than above for orbits with eccentricities greater than 0.3. This // than above for orbits with eccentricities greater than 0.3. This
// is from Jean Meeus's _Astronomical Algorithms_ (2nd ed), p. 199 // is from Jean Meeus's _Astronomical Algorithms_ (2nd ed), p. 199
struct SolveKeplerFunc2 : public unary_function<double, double> struct SolveKeplerFunc2 : public unary_function<double, double>
{ {
double ecc; double ecc;
double M; double M;
SolveKeplerFunc2(double _ecc, double _M) : ecc(_ecc), M(_M) {} SolveKeplerFunc2(double _ecc, double _M) : ecc(_ecc), M(_M) {}
double operator()(double x) const double operator()(double x) const
{ {
return x + (M + ecc * sin(x) - x) / (1 - ecc * cos(x)); return x + (M + ecc * sin(x) - x) / (1 - ecc * cos(x));
} }
}; };
double sign(double x) double sign(double x)
{ {
if (x < 0.) if (x < 0.)
return -1.; return -1.;
else if (x > 0.) else if (x > 0.)
return 1.; return 1.;
else else
return 0.; return 0.;
} }
struct SolveKeplerLaguerreConway : public unary_function<double, double> struct SolveKeplerLaguerreConway : public unary_function<double, double>
{ {
double ecc; double ecc;
double M; double M;
SolveKeplerLaguerreConway(double _ecc, double _M) : ecc(_ecc), M(_M) {} SolveKeplerLaguerreConway(double _ecc, double _M) : ecc(_ecc), M(_M)
// cf Heafner, Fundamental Ephemeris Computations p.73 {}
// GZ: note&add Heafner's initial guess for E! // cf Heafner, Fundamental Ephemeris Computations p.73
double operator()(double E) const // GZ: note&add Heafner's initial guess for E!
{ double operator()(double E) const
double s = ecc * sin(E); {
double c = ecc * cos(E); double s = ecc * sin(E);
double f = E - s - M; double c = ecc * cos(E);
double f1 = 1 - c; double f = E - s - M;
double f2 = s; double f1 = 1 - c;
E += -5 * f / (f1 + sign(f1) * sqrt(abs(16 * f1 * f1 - 20 * f * f2) double f2 = s;
)); E += -5 * f / (f1 + sign(f1) * std::sqrt(abs(16 * f1 * f1 -
20 * f * f2)));
return E; return E;
} }
}; };
struct SolveKeplerLaguerreConwayHyp : public unary_function<double, double> struct SolveKeplerLaguerreConwayHyp : public unary_function<double, double>
{ {
double ecc; double ecc;
double M; double M;
SolveKeplerLaguerreConwayHyp(double _ecc, double _M) : ecc(_ecc), M(_M) SolveKeplerLaguerreConwayHyp(double _ecc, double _M) : ecc(_ecc), M(
{} _M) {}
// cf Heafner, Fundamental Ephemeris Computations p.73 // cf Heafner, Fundamental Ephemeris Computations p.73
double operator()(double x) const double operator()(double x) const
{ {
double s = ecc * sinh(x); double s = ecc * sinh(x);
double c = ecc * cosh(x); double c = ecc * cosh(x);
double f = s - x - M; double f = s - x - M;
double f1 = c - 1; double f1 = c - 1;
double f2 = s; double f2 = s;
x += -5 * f / (f1 + sign(f1) * sqrt(abs(16 * f1 * f1 - 20 * f * f2) x += -5 * f / (f1 + sign(f1) * std::sqrt(abs(16 * f1 * f1 -
)); 20 * f * f2)));
return x; return x;
} }
}; };
typedef pair<double, double> Solution; typedef pair<double, double> Solution;
double EllipticalOrbit::eccentricAnomaly(const double M) const double EllipticalOrbit::eccentricAnomaly(const double M) const
{ {
if (eccentricity == 0.0) if (eccentricity == 0.0)
{ {
// Circular orbit // Circular orbit
return M; return M;
} }
else if (eccentricity < 0.2) else if (eccentricity < 0.2)
{ {
// Low eccentricity, so use the standard iteration technique // Low eccentricity, so use the standard iteration technique
Solution sol = solveIteration_fixed(SolveKeplerFunc1(eccentricity, Solution sol = solveIteration_fixed(SolveKeplerFunc1(eccentr
M), M, 5); icity, M), M, 5);
return sol.first; return sol.first;
} }
else if (eccentricity < 0.9) else if (eccentricity < 0.9)
{ {
// Higher eccentricity elliptical orbit; use a more complex but // Higher eccentricity elliptical orbit; use a more complex
// much faster converging iteration. but
Solution sol = solveIteration_fixed(SolveKeplerFunc2(eccentricity, // much faster converging iteration.
M), M, 6); Solution sol = solveIteration_fixed(SolveKeplerFunc2(eccentr
// Debugging icity, M), M, 6);
// qDebug("ecc: %f, error: %f mas\n", // Debugging
// eccentricity, radToDeg(sol.second) * 3600000); // qDebug("ecc: %f, error: %f mas\n",
return sol.first; // eccentricity, radToDeg(sol.second) * 3600000);
} return sol.first;
else if (eccentricity < 1.0) }
{ else if (eccentricity < 1.0)
// Extremely stable Laguerre-Conway method for solving Kepler's {
// equation. Only use this for high-eccentricity orbits, as it // Extremely stable Laguerre-Conway method for solving Keple
// requires more calculation. r's
double E = M + 0.85 * eccentricity * sign(sin(M)); // equation. Only use this for high-eccentricity orbits, as
Solution sol = solveIteration_fixed(SolveKeplerLaguerreConway(eccen it
tricity, M), E, 8); // requires more calculation.
return sol.first; double E = M + 0.85 * eccentricity * sign(sin(M));
} Solution sol = solveIteration_fixed(SolveKeplerLaguerreConwa
else if (eccentricity == 1.0) y(eccentricity, M), E, 8);
{ return sol.first;
// parabolic orbit; very common for comets }
// TODO: handle this. else if (eccentricity == 1.0)
// Problem: E does not make sense here. True anomaly quantities (rS {
inNu, rCosNu) computed directly. // parabolic orbit; very common for comets
// Anyhow, Comets use CometOrbit class. // TODO: handle this.
return M; // Problem: E does not make sense here. True anomaly quantit
} ies (rSinNu, rCosNu) computed directly.
else // Anyhow, Comets use CometOrbit class.
{ return M;
// Laguerre-Conway method for hyperbolic (ecc > 1) orbits. }
double E = log(2 * M / eccentricity + 1.85); else
Solution sol = solveIteration_fixed(SolveKeplerLaguerreConwayHyp(ec {
centricity, M), E, 30); // Laguerre-Conway method for hyperbolic (ecc > 1) orbits.
return sol.first; double E = log(2 * M / eccentricity + 1.85);
} Solution sol = solveIteration_fixed(SolveKeplerLaguerreConwa
yHyp(eccentricity, M), E, 30);
return sol.first;
}
} }
Vec3d EllipticalOrbit::positionAtE(const double E) const Vec3d EllipticalOrbit::positionAtE(const double E) const
{ {
double x, z; double x, z;
if (eccentricity < 1.0) if (eccentricity < 1.0)
{ {
double a = pericenterDistance / (1.0 - eccentricity); double a = pericenterDistance / (1.0 - eccentricity);
x = a * (cos(E) - eccentricity); x = a * (cos(E) - eccentricity);
z = a * sqrt(1 - eccentricity * eccentricity) * -sin(E); z = a * std::sqrt(1 - eccentricity * eccentricity) * -sin(E)
} ;
else if (eccentricity > 1.0) // N.B. This is odd at least: elliptical m }
ust have ecc<1! else if (eccentricity > 1.0) // N.B. This is odd at least: elliptica
{ l must have ecc<1!
double a = pericenterDistance / (1.0 - eccentricity); {
x = -a * (eccentricity - cosh(E)); double a = pericenterDistance / (1.0 - eccentricity);
z = -a * sqrt(eccentricity * eccentricity - 1) * -sinh(E); x = -a * (eccentricity - cosh(E));
} z = -a * std::sqrt(eccentricity * eccentricity - 1) * -sinh(
else E);
{ }
// TODO: Handle parabolic orbits else
x = 0.0; {
z = 0.0; // TODO: Handle parabolic orbits
} x = 0.0;
z = 0.0;
Mat4d R = (Mat4d::zrotation(ascendingNode) * }
Mat4d::xrotation(inclination) *
Mat4d::zrotation(argOfPeriapsis)); Mat4d R = (Mat4d::zrotation(ascendingNode) *
Mat4d::xrotation(inclination) *
Mat4d::zrotation(argOfPeriapsis));
return R * Vec3d(x, -z, 0); return R * Vec3d(x, -z, 0);
} }
// Return the offset from the center. // Return the offset from the center.
Vec3d EllipticalOrbit::positionAtTime(const double JD) const Vec3d EllipticalOrbit::positionAtTime(const double JD) const
{ {
double meanMotion = 2.0 * M_PI / period; double meanMotion = 2.0 * M_PI / period;
double meanAnomaly = meanAnomalyAtEpoch + (JD-epoch) * meanMotion; double meanAnomaly = meanAnomalyAtEpoch + (JD-epoch) * meanMotion;
double E = eccentricAnomaly(meanAnomaly); double E = eccentricAnomaly(meanAnomaly);
return positionAtE(E); return positionAtE(E);
} }
//void EllipticalOrbit::positionAtTime(double JD, double * X, double * Y, d ouble * Z) const //void EllipticalOrbit::positionAtTime(double JD, double * X, double * Y, d ouble * Z) const
//{ //{
// Vec3d pos = positionAtTime(JD); // Vec3d pos = positionAtTime(JD);
// *X=pos[2]; // *X=pos[2];
// *Y=pos[0]; // *Y=pos[0];
// *Z=pos[1]; // *Z=pos[1];
//} //}
//void EllipticalOrbit::positionAtTimev(double JD, double* v) //void EllipticalOrbit::positionAtTimev(double JD, double* v)
//{ //{
// Vec3d pos = positionAtTime(JD); // Vec3d pos = positionAtTime(JD);
// v[0]=pos[2]; // v[0]=pos[2];
// v[1]=pos[0]; // v[1]=pos[0];
// v[2]=pos[1]; // v[2]=pos[1];
//} //}
void EllipticalOrbit::positionAtTimevInVSOP87Coordinates(const double JD, d ouble* v) const void EllipticalOrbit::positionAtTimevInVSOP87Coordinates(const double JD, d ouble* v) const
{ {
Vec3d pos = positionAtTime(JD); Vec3d pos = positionAtTime(JD);
v[0] = rotateToVsop87[0]*pos[0] + rotateToVsop87[1]*pos[1] + rotateToVsop v[0] = rotateToVsop87[0]*pos[0] + rotateToVsop87[1]*pos[1] + rotateT
87[2]*pos[2]; oVsop87[2]*pos[2];
v[1] = rotateToVsop87[3]*pos[0] + rotateToVsop87[4]*pos[1] + rotateToVsop v[1] = rotateToVsop87[3]*pos[0] + rotateToVsop87[4]*pos[1] + rotateT
87[5]*pos[2]; oVsop87[5]*pos[2];
v[2] = rotateToVsop87[6]*pos[0] + rotateToVsop87[7]*pos[1] + rotateToVsop v[2] = rotateToVsop87[6]*pos[0] + rotateToVsop87[7]*pos[1] + rotateT
87[8]*pos[2]; oVsop87[8]*pos[2];
} }
double EllipticalOrbit::getPeriod() const double EllipticalOrbit::getPeriod() const
{ {
return period; return period;
} }
double EllipticalOrbit::getBoundingRadius() const double EllipticalOrbit::getBoundingRadius() const
{ {
// TODO: watch out for unbounded parabolic and hyperbolic orbits // TODO: watch out for unbounded parabolic and hyperbolic orbits
return pericenterDistance * ((1.0 + eccentricity) / (1.0 - eccentricity return pericenterDistance * ((1.0 + eccentricity) / (1.0 - eccentric
)); ity));
} }
void EllipticalOrbit::sample(double, double, int nSamples, OrbitSampleProc& proc) const void EllipticalOrbit::sample(double, double, int nSamples, OrbitSampleProc& proc) const
{ {
double dE = 2. * M_PI / (double) nSamples; double dE = 2. * M_PI / (double) nSamples;
for (int i = 0; i < nSamples; i++) for (int i = 0; i < nSamples; i++)
proc.sample(positionAtE(dE * i)); proc.sample(positionAtE(dE * i));
} }
Vec3d CachingOrbit::positionAtTime(double jd) const Vec3d CachingOrbit::positionAtTime(double jd) const
{ {
if (jd != lastTime) if (jd != lastTime)
{ {
lastTime = jd; lastTime = jd;
lastPosition = computePosition(jd); lastPosition = computePosition(jd);
} }
return lastPosition; return lastPosition;
} }
void CachingOrbit::sample(double start, double t, int nSamples, void CachingOrbit::sample(double start, double t, int nSamples, OrbitSample
OrbitSampleProc& proc) const Proc& proc) const
{ {
double dt = t / (double) nSamples; double dt = t / (double) nSamples;
for (int i = 0; i < nSamples; i++) for (int i = 0; i < nSamples; i++)
proc.sample(positionAtTime(start + dt * i)); proc.sample(positionAtTime(start + dt * i));
} }
 End of changes. 32 change blocks. 
359 lines changed or deleted 381 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/