Orbit.cpp   Orbit.cpp
skipping to change at line 130 skipping to change at line 130
double M = fmod(n*dt,2*M_PI); // Mean Anomaly double M = fmod(n*dt,2*M_PI); // Mean Anomaly
if (M < 0.0) M += 2.0*M_PI; if (M < 0.0) M += 2.0*M_PI;
// double E = M; // double E = M;
// for (;;) { // Newton(?) Solve Kepler's equation (similar to Meeus se cond method, Astro.Alg. 1998 p.199) // for (;;) { // Newton(?) Solve Kepler's equation (similar to Meeus se cond method, Astro.Alg. 1998 p.199)
// const double Ep = E; // const double Ep = E;
// E -= (M-E+e*sin(E))/(e*cos(E)-1); // E -= (M-E+e*sin(E))/(e*cos(E)-1);
// if (fabs(E-Ep) < EPSILON) break; // if (fabs(E-Ep) < EPSILON) break;
// } // }
// GZ: Comet orbits are quite often near-parabolic, where this may stil l only converge slowly. // GZ: Comet orbits are quite often near-parabolic, where this may stil l only converge slowly.
// Better always use Laguerre-Conway. See Heafner, Ch. 5.3 // Better always use Laguerre-Conway. See Heafner, Ch. 5.3
// Ouch! https://bugs.launchpad.net/stellarium/+bug/1465112 ==>It seems
we still need an escape counter!
// Debug line in test case fabs(E-Ep) indicates it usually takes 2-3,
occasionally up to 6 cycles.
// It seems safe to assume 10 should not be exceeded. N.B.: A GPU fixed
-loopcount implementation could go for 8 passes.
double E=M+0.85*e*StelUtils::sign(sin(M)); double E=M+0.85*e*StelUtils::sign(sin(M));
int escape=0;
for (;;) for (;;)
{ {
const double Ep=E; const double Ep=E;
const double f2=e*sin(E); const double f2=e*sin(E);
const double f=E-f2-M; const double f=E-f2-M;
const double f1=1.0-e*cos(E); const double f1=1.0-e*cos(E);
E+= (-5.0*f)/(f1+StelUtils::sign(f1)*std::sqrt(fabs(16.0*f1* f1-20.0*f*f2))); E+= (-5.0*f)/(f1+StelUtils::sign(f1)*std::sqrt(fabs(16.0*f1* f1-20.0*f*f2)));
if (fabs(E-Ep) < EPSILON) break; if (fabs(E-Ep) < EPSILON)
{
//qDebug() << "Ell. orbit with eccentricity " << e <
< "Escaping after" << escape << "loops at E-Ep=" << E-Ep;
break;
}
escape++;
if (escape>10)
{
qDebug() << "Ell. orbit with eccentricity " << e <<
"would have caused endless loop. Escaping after 10 runs at E-Ep=" << E-Ep;
break;
}
} }
// Note: q=a*(1-e) // Note: q=a*(1-e)
const double h1 = q*std::sqrt((1.0+e)/(1.0-e)); // elsewhere: a sqr t(1-e²) ... q / (1-e) sqrt( (1+e)(1-e)) = q sqrt((1+e)/(1-e)) const double h1 = q*std::sqrt((1.0+e)/(1.0-e)); // elsewhere: a sqr t(1-e²) ... q / (1-e) sqrt( (1+e)(1-e)) = q sqrt((1+e)/(1-e))
rCosNu = a*(cos(E)-e); rCosNu = a*(cos(E)-e);
rSinNu = h1*sin(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
skipping to change at line 223 skipping to change at line 237
orbitGood(orbitGoodDays) orbitGood(orbitGoodDays)
{ {
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;
rotateToVsop87[2] = s_nod*s_obl; rotateToVsop87[2] = s_nod*s_obl;
rotateToVsop87[3] = s_nod*cj+c_nod*c_obl*sj; rotateToVsop87[3] = s_nod*cj+c_nod*c_obl*sj;
rotateToVsop87[4] = -s_nod*sj+c_nod*c_obl*cj; rotateToVsop87[4] = -s_nod*sj+c_nod*c_obl*cj;
rotateToVsop87[5] = -c_nod*s_obl; rotateToVsop87[5] = -c_nod*s_obl;
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 JDool updateVelocityVector) void CometOrbit::positionAtTimevInVSOP87Coordinates(double JDE, double *v, bool updateVelocityVector)
{ {
JD -= t0; JDE -= t0;
double rCosNu,rSinNu; double rCosNu,rSinNu;
// temporary solve freezes for near-parabolic comets - using (e < 0.999 9) for elliptical orbits // temporary solve freezes for near-parabolic comets - using (e < 0.999 9) 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,JDith Laguerre-Conway, I dare to go for 1.0. if (e < 1.0) InitEll(q,n,e,JDE,rCosNu,rSinNu); // GZ: After solving with Laguerre-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,JDE,rCosNu,rSinNu);
} }
else InitPar(q,n,JD,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;
v[2] = rotateToVsop87[6]*p0 + rotateToVsop87[7]*p1 + rotateToVsop87[ 8]*p2; v[2] = rotateToVsop87[6]*p0 + rotateToVsop87[7]*p1 + rotateToVsop87[ 8]*p2;
if (updateVelocityVector) if (updateVelocityVector)
{ {
rdot.set(s0, s1, s2); rdot.set(s0, s1, s2);
updateTails=true; updateTails=true;
skipping to change at line 472 skipping to change at line 477
} }
Mat4d R = (Mat4d::zrotation(ascendingNode) * Mat4d R = (Mat4d::zrotation(ascendingNode) *
Mat4d::xrotation(inclination) * Mat4d::xrotation(inclination) *
Mat4d::zrotation(argOfPeriapsis)); 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 JDE) const
{ {
double meanMotion = 2.0 * M_PI / period; double meanMotion = 2.0 * M_PI / period;
double meanAnomaly = meanAnomalyAtEpoch + (JD-epoch) * meanMotion; double meanAnomaly = meanAnomalyAtEpoch + (JDE-epoch) * meanMotion;
double E = eccentricAnomaly(meanAnomaly); double E = eccentricAnomaly(meanAnomaly);
return positionAtE(E); return positionAtE(E);
} }
//void EllipticalOrbit::positionAtTime(double JDouble * Z) const //void EllipticalOrbit::positionAtTime(double JDE, double * X, double * Y, double * Z) const
//{ //{
// Vec3d pos = positionAtTime(JD); // Vec3d pos = positionAtTime(JDE);
// *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 JDE, double* v)
//{ //{
// Vec3d pos = positionAtTime(JD); // Vec3d pos = positionAtTime(JDE);
// 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 JDouble* v) const void EllipticalOrbit::positionAtTimevInVSOP87Coordinates(const double JDE, double* v) const
{ {
Vec3d pos = positionAtTime(JD); Vec3d pos = positionAtTime(JDE);
v[0] = rotateToVsop87[0]*pos[0] + rotateToVsop87[1]*pos[1] + rotateT oVsop87[2]*pos[2]; v[0] = rotateToVsop87[0]*pos[0] + rotateToVsop87[1]*pos[1] + rotateT oVsop87[2]*pos[2];
v[1] = rotateToVsop87[3]*pos[0] + rotateToVsop87[4]*pos[1] + rotateT oVsop87[5]*pos[2]; v[1] = rotateToVsop87[3]*pos[0] + rotateToVsop87[4]*pos[1] + rotateT oVsop87[5]*pos[2];
v[2] = rotateToVsop87[6]*pos[0] + rotateToVsop87[7]*pos[1] + rotateT oVsop87[8]*pos[2]; v[2] = rotateToVsop87[6]*pos[0] + rotateToVsop87[7]*pos[1] + rotateT oVsop87[8]*pos[2];
} }
double EllipticalOrbit::getPeriod() const double EllipticalOrbit::getPeriod() const
{ {
return period; return period;
} }
skipping to change at line 523 skipping to change at line 528
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));
} }
Vec3d CachingOrbit::positionAtTime(double ) const Vec3d CachingOrbit::positionAtTime(double JDE) const
{ {
if ( != lastTime) if (JDE != lastTime)
{ {
lastTime = lastTime = JDE;
lastPosition = 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. 20 change blocks.
27 lines changed or deleted 37 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/