Orbit.cpp   Orbit.cpp 
skipping to change at line 157 skipping to change at line 157
//! @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 &rdotz, const bool withVelVector=false, const double e=0.0, const d ouble q=0.0) { double &rx,double &ry,double &rz, double &rdotx, double &rdoty, double &rdotz, const bool withVelVector=false, const double e=0.0, const d ouble q=0.0) {
// qDebug() << "Init3D";
const double cw = cos(w); const double cw = cos(w);
const double sw = sin(w); const double sw = sin(w);
const double cOm = cos(Omega); const double cOm = cos(Omega);
const double sOm = sin(Omega); const double sOm = sin(Omega);
const double ci = cos(i); const double ci = cos(i);
const double si = sin(i); const double si = sin(i);
const double Px=-sw*sOm*ci+cw*cOm; // Heafner, 5.3.1 Px const double Px=-sw*sOm*ci+cw*cOm; // Heafner, 5.3.1 Px
const double Qx=-cw*sOm*ci-sw*cOm; // Heafner, 5.3.4 Qx const double Qx=-cw*sOm*ci-sw*cOm; // Heafner, 5.3.4 Qx
const double Py= sw*cOm*ci+cw*sOm; // Heafner, 5.3.2 Py const double Py= sw*cOm*ci+cw*sOm; // Heafner, 5.3.2 Py
const double Qy= cw*cOm*ci-sw*sOm; // Heafner, 5.3.5 Qy const double Qy= cw*cOm*ci-sw*sOm; // Heafner, 5.3.5 Qy
skipping to change at line 193 skipping to change at line 192
} }
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, // GZ: I don't see a double parentRotObliquity,
ny use for this, should be 0 for comets. double parentRotAscendingnode,
double parentRotAscendingnode, // GZ: I don't see a double parentRotJ2000Longitude)
ny use for this, should be 0 for comets.
double parentRotJ2000Longitude) // GZ: I don't see a
ny use for this, should be 0 for comets.
//) //)
:q(pericenterDistance),e(eccentricity),i(inclination), :q(pericenterDistance),e(eccentricity),i(inclination),
Om(ascendingNode),w(argOfPerhelion),t0(timeAtPerihelion), Om(ascendingNode),w(argOfPerhelion),t0(timeAtPerihelion),
n(meanMotion), updateTails(true), orbitGood(orbitGoodDays) { n(meanMotion), updateTails(true), orbitGood(orbitGoodDays) {
// qDebug() << "CometOrbit::()";
rdot.set(0.0, 0.0, 0.0); rdot.set(0.0, 0.0, 0.0);
const double c_obl = cos(parentRotObliquity); // 1? const double c_obl = cos(parentRotObliquity);
const double s_obl = sin(parentRotObliquity); // 0? const double s_obl = sin(parentRotObliquity);
const double c_nod = cos(parentRotAscendingnode); // 1? const double c_nod = cos(parentRotAscendingnode);
const double s_nod = sin(parentRotAscendingnode); // 0? const double s_nod = sin(parentRotAscendingnode);
const double cj = cos(parentRotJ2000Longitude); // 1? const double cj = cos(parentRotJ2000Longitude);
const double sj = sin(parentRotJ2000Longitude); // 0? const double sj = sin(parentRotJ2000Longitude);
// GZ: Test my assumptions before breaking anything...
// Q_ASSERT(c_obl==1.0);
// Q_ASSERT(c_nod==1.0);
// Q_ASSERT(cj ==1.0);
// Q_ASSERT(s_obl==0.0);
// Q_ASSERT(s_nod==0.0);
// Q_ASSERT(sj ==0.0);
// GZ NO, this is necessary!
// rotateToVsop87[0] = c_nod; // rotateToVsop87[0] = c_nod;
// rotateToVsop87[1] = -s_nod * c_obl; // rotateToVsop87[1] = -s_nod * c_obl;
// rotateToVsop87[2] = s_nod * s_obl; // rotateToVsop87[2] = s_nod * s_obl;
// rotateToVsop87[3] = s_nod; // rotateToVsop87[3] = s_nod;
// rotateToVsop87[4] = c_nod * c_obl; // rotateToVsop87[4] = c_nod * c_obl;
// rotateToVsop87[5] = -c_nod * s_obl; // rotateToVsop87[5] = -c_nod * s_obl;
// rotateToVsop87[6] = 0.0; // rotateToVsop87[6] = 0.0;
// rotateToVsop87[7] = s_obl; // rotateToVsop87[7] = s_obl;
// rotateToVsop87[8] = c_obl; // rotateToVsop87[8] = c_obl;
rotateToVsop87[0] = c_nod*cj-s_nod*c_obl*sj; // 1 in case of comet orbit rotateToVsop87[0] = c_nod*cj-s_nod*c_obl*sj;
ing the sun, these names are misleading, however they do the right thing. rotateToVsop87[1] = -c_nod*sj-s_nod*c_obl*cj;
rotateToVsop87[1] = -c_nod*sj-s_nod*c_obl*cj; // 0 OK, this is NOT an ide rotateToVsop87[2] = s_nod*s_obl;
ntity matrix... rotateToVsop87[3] = s_nod*cj+c_nod*c_obl*sj;
rotateToVsop87[2] = s_nod*s_obl; // 0 rotateToVsop87[4] = -s_nod*sj+c_nod*c_obl*cj;
rotateToVsop87[3] = s_nod*cj+c_nod*c_obl*sj; // 0 rotateToVsop87[5] = -c_nod*s_obl;
rotateToVsop87[4] = -s_nod*sj+c_nod*c_obl*cj; // 1 rotateToVsop87[6] = s_obl*sj;
rotateToVsop87[5] = -c_nod*s_obl; // 0 rotateToVsop87[7] = s_obl*cj;
rotateToVsop87[6] = s_obl*sj; // 0 rotateToVsop87[8] = c_obl;
rotateToVsop87[7] = s_obl*cj; // 0
rotateToVsop87[8] = c_obl; // 1
// qDebug() << "CometOrbit::()...done"; // qDebug() << "CometOrbit::()...done";
} }
void CometOrbit::positionAtTimevInVSOP87Coordinates(double JD, double *v, b ool updateVelocityVector) { void CometOrbit::positionAtTimevInVSOP87Coordinates(double JD, double *v, b ool updateVelocityVector) {
JD -= t0; JD -= t0;
double rCosNu,rSinNu; double rCosNu,rSinNu;
// temporary solve freezes for near-parabolic comets - using (e < 0.9999) for elliptical orbits // temporary solve freezes for near-parabolic comets - using (e < 0.9999) for elliptical orbits
// TODO: improve calculations orbits for near-parabolic comets --AW // TODO: improve calculations orbits for near-parabolic comets --AW
// if (e < 0.9999) InitEll(q,n,e,JD,a1,a2); // if (e < 0.9999) InitEll(q,n,e,JD,a1,a2);
if (e < 1.0) InitEll(q,n,e,JD,rCosNu,rSinNu); // GZ: After solving with L aguerre-Conway, I dare to go for 1.0. if (e < 1.0) InitEll(q,n,e,JD,rCosNu,rSinNu); // GZ: After solving with L aguerre-Conway, I dare to go for 1.0.
else if (e > 1.0) { else if (e > 1.0) {
// qDebug() << "Hyperbolic orbit for ecc=" << e << ", i=" << i << ", w =" << w << ", Mean Motion n=" << n; // qDebug() << "Hyperbolic orbit for ecc=" << e << ", i=" << i << ", w =" << w << ", Mean Motion n=" << n;
InitHyp(q,n,e,JD,rCosNu,rSinNu); InitHyp(q,n,e,JD,rCosNu,rSinNu);
} }
else InitPar(q,n,JD,rCosNu,rSinNu); else InitPar(q,n,JD,rCosNu,rSinNu);
double p0,p1,p2, s0, s1, s2; double p0,p1,p2, s0, s1, s2;
Init3D(i,Om,w,rCosNu,rSinNu,p0,p1,p2, s0, s1, s2, updateVelocityVector, e , q); Init3D(i,Om,w,rCosNu,rSinNu,p0,p1,p2, s0, s1, s2, updateVelocityVector, e , q);
// GZ: The next 3 lines are meaningless for a comet orbiting the sun. Or not?
v[0] = rotateToVsop87[0]*p0 + rotateToVsop87[1]*p1 + rotateToVsop87[2]*p2 ; v[0] = rotateToVsop87[0]*p0 + rotateToVsop87[1]*p1 + rotateToVsop87[2]*p2 ;
v[1] = rotateToVsop87[3]*p0 + rotateToVsop87[4]*p1 + rotateToVsop87[5]*p2 ; v[1] = rotateToVsop87[3]*p0 + rotateToVsop87[4]*p1 + rotateToVsop87[5]*p2 ;
v[2] = rotateToVsop87[6]*p0 + rotateToVsop87[7]*p1 + rotateToVsop87[8]*p2 ; v[2] = rotateToVsop87[6]*p0 + rotateToVsop87[7]*p1 + rotateToVsop87[8]*p2 ;
//GZ replace with:
//v[0]=p0;
//v[1]=p1;
//v[2]=p2;
if (updateVelocityVector) { if (updateVelocityVector) {
rdot.set(s0, s1, s2); // TODO: ROTATE? rdot.set(s0, s1, s2);
updateTails=true; 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,
skipping to change at line 311 skipping to change at line 298
rotateToVsop87[7] = s_obl*cj; rotateToVsop87[7] = s_obl*cj;
rotateToVsop87[8] = c_obl; 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.)
skipping to change at line 350 skipping to change at line 337
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 // cf Heafner, Fundamental Ephemeris Computations p.73
// GZ: note&add Heafner's initial guess for E! // GZ: note&add Heafner's initial guess for E!
double operator()(double E) const double operator()(double E) const
{ {
double s = ecc * sin(E); double s = ecc * sin(E);
double c = ecc * cos(E); double c = ecc * cos(E);
double f = E - s - M; double f = E - s - M;
double f1 = 1 - c; double f1 = 1 - c;
double f2 = s; double f2 = s;
E += -5 * f / (f1 + sign(f1) * sqrt(abs(16 * f1 * f1 - 20 * f * f2) )); E += -5 * f / (f1 + sign(f1) * 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) * sqrt(abs(16 * f1 * f1 - 20 * f * f2) ));
skipping to change at line 415 skipping to change at line 402
Solution sol = solveIteration_fixed(SolveKeplerFunc2(eccentricity, M), M, 6); Solution sol = solveIteration_fixed(SolveKeplerFunc2(eccentricity, M), M, 6);
// Debugging // Debugging
// qDebug("ecc: %f, error: %f mas\n", // qDebug("ecc: %f, error: %f mas\n",
// eccentricity, radToDeg(sol.second) * 3600000); // eccentricity, radToDeg(sol.second) * 3600000);
return sol.first; return sol.first;
} }
else if (eccentricity < 1.0) else if (eccentricity < 1.0)
{ {
// Extremely stable Laguerre-Conway method for solving Kepler's // Extremely stable Laguerre-Conway method for solving Kepler's
// equation. Only use this for high-eccentricity orbits, as it // equation. Only use this for high-eccentricity orbits, as it
// requires more calcuation. // requires more calculation.
double E = M + 0.85 * eccentricity * sign(sin(M)); double E = M + 0.85 * eccentricity * sign(sin(M));
Solution sol = solveIteration_fixed(SolveKeplerLaguerreConway(eccen tricity, M), E, 8); Solution sol = solveIteration_fixed(SolveKeplerLaguerreConway(eccen tricity, M), E, 8);
return sol.first; return sol.first;
} }
else if (eccentricity == 1.0) else if (eccentricity == 1.0)
{ {
// parabolic orbit; very common for comets // parabolic orbit; very common for comets
// TODO: handle this. // TODO: handle this.
// Problem: E does not make sense here. True anomaly quantities (rS inNu, rCosNu) computed directly. // Problem: E does not make sense here. True anomaly quantities (rS inNu, rCosNu) computed directly.
// Anyhow, Comets use CometOrbit class. // Anyhow, Comets use CometOrbit class.
 End of changes. 13 change blocks. 
44 lines changed or deleted 26 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/