Orbit.cpp   Orbit.cpp 
skipping to change at line 102 skipping to change at line 102
a2 = 2.0*q*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 &rCosNu, double &rSinNu) static void InitPar(const double q, const double n, const double dt, double &rCosNu, double &rSinNu)
{ {
// qDebug() << "InitPar"; // qDebug() << "InitPar";
// const double M=dt*sqrt(GAUSS_GRAV_CONST/(2.0*q*q*q)); // const double M=dt*sqrt(GAUSS_GRAV_CONST/(2.0*q*q*q));
// const double W=1.5*M; // const double W=1.5*M;
const double W=dt*n; const double W=dt*n;
const double Y=cbrt(W+std::sqrt(W*W+1)); const double Y=cbrt(W+std::sqrt(W*W+1.));
const double tanNu2=Y-1.0/Y; // Heafner (5.5.8) has an error here, w rites (Y-1)/Y. const double tanNu2=Y-1.0/Y; // Heafner (5.5.8) has an error here, w rites (Y-1)/Y.
rCosNu=q*(1.0-tanNu2*tanNu2); rCosNu=q*(1.0-tanNu2*tanNu2);
rSinNu=2.0*q*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
skipping to change at line 218 skipping to change at line 218
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), : q(pericenterDistance),
e(eccentricity), e(eccentricity),
i(inclination), i(inclination),
Om(ascendingNode), Om(ascendingNode),
w(argOfPerhelion), w(argOfPerhelion),
t0(timeAtPerihelion), t0(timeAtPerihelion),
n(meanMotion), n(meanMotion),
updateTails(true), updateTails(true),
orbitGood(orbitGoodDays) orbitGood(orbitGoodDays)
{ {
// GZ MAKE SURE THIS IS ALWAYS 0/0/0. ==> OK.
//qDebug() << "parentRotObliquity" << parentRotObliquity << "parentR
otAscendingnode" << parentRotAscendingnode << "parentRotJ2000Longitude" <<
parentRotJ2000Longitude;
rdot.set(0.0, 0.0, 0.0); rdot.set(0.0, 0.0, 0.0);
const double c_obl = cos(parentRotObliquity); const double c_obl = cos(parentRotObliquity);
const double s_obl = sin(parentRotObliquity); const double s_obl = sin(parentRotObliquity);
const double c_nod = cos(parentRotAscendingnode); const double c_nod = cos(parentRotAscendingnode);
const double s_nod = sin(parentRotAscendingnode); const double s_nod = sin(parentRotAscendingnode);
const double cj = cos(parentRotJ2000Longitude); const double cj = cos(parentRotJ2000Longitude);
const double sj = sin(parentRotJ2000Longitude); const double sj = sin(parentRotJ2000Longitude);
rotateToVsop87[0] = c_nod*cj-s_nod*c_obl*sj; rotateToVsop87[0] = c_nod*cj-s_nod*c_obl*sj;
rotateToVsop87[1] = -c_nod*sj-s_nod*c_obl*cj; rotateToVsop87[1] = -c_nod*sj-s_nod*c_obl*cj;
skipping to change at line 253 skipping to change at line 254
rotateToVsop87[6] = s_obl*sj; rotateToVsop87[6] = s_obl*sj;
rotateToVsop87[7] = s_obl*cj; rotateToVsop87[7] = s_obl*cj;
rotateToVsop87[8] = c_obl; rotateToVsop87[8] = c_obl;
// qDebug() << "CometOrbit::()...done"; // qDebug() << "CometOrbit::()...done";
} }
void CometOrbit::positionAtTimevInVSOP87Coordinates(double JDE, double *v, bool updateVelocityVector) void CometOrbit::positionAtTimevInVSOP87Coordinates(double JDE, double *v, bool updateVelocityVector)
{ {
JDE -= t0; JDE -= t0;
double rCosNu,rSinNu; double rCosNu,rSinNu;
if (e < 1.0) InitEll(q,n,e,JDE,rCosNu,rSinNu); // GZ: After solving with Laguerre-Conway, I dare to go for 1.0. if (e < 1.0) InitEll(q,n,e,JDE,rCosNu,rSinNu); // Laguerre-Conway se ems stable enough 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,JDE,rCosNu,rSinNu); InitHyp(q,n,e,JDE,rCosNu,rSinNu);
} }
else InitPar(q,n,JDE,rCosNu,rSinNu); else InitPar(q,n,JDE,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, updateVelocityVect or, e, q); Init3D(i,Om,w,rCosNu,rSinNu,p0,p1,p2, s0, s1, s2, updateVelocityVect or, e, q);
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;
skipping to change at line 525 skipping to change at line 526
return pericenterDistance * ((1.0 + eccentricity) / (1.0 - eccentric ity)); 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));
} }
/*
* Stuff found unused and deactivated pre-0.15
Vec3d CachingOrbit::positionAtTime(double JDE) const Vec3d CachingOrbit::positionAtTime(double JDE) const
{ {
if (JDE != lastTime) if (JDE != lastTime)
{ {
lastTime = JDE; lastTime = JDE;
lastPosition = computePosition(JDE); lastPosition = computePosition(JDE);
} }
return lastPosition; return lastPosition;
} }
void CachingOrbit::sample(double start, double t, int nSamples, OrbitSample Proc& proc) const void CachingOrbit::sample(double start, double t, int nSamples, OrbitSample 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. 6 change blocks. 
3 lines changed or deleted 8 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/