Orbit.cpp   Orbit.cpp 
// orbit.cpp // orbit.cpp
// //
// Copyright (C) 2001, Chris Laurel <claurel@shatters.net> // Copyright (C) 2001, Chris Laurel <claurel@shatters.net>
// //
// CometOrbit, InitHyp,InitPar,InitEll,Init3D: Copyright (C) 1995,2007,2008 Johannes Gajdosik // CometOrbit, InitHyp,InitPar,InitEll,Init3D: Copyright (C) 1995,2007,2008 Johannes Gajdosik
// InitHyp,InitPar,InitEll,Init3D checked, improved, extended, annotated 20
13 Georg Zotti (GZ).
// Algorithms identified and extended from:
// Meeus: Astronomical Algorithms 1998
// Heafner: Fundamental Ephemeris Computations 1999
// //
// This program is free software; you can redistribute it and/or // This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License // modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 2 // as published by the Free Software Foundation; either version 2
// of the License, or (at your option) any later version. // of the License, or (at your option) any later version.
#include "Solve.hpp"
#include "Orbit.hpp"
#include "StelUtils.hpp"
#include <functional> #include <functional>
#include <algorithm> #include <algorithm>
#include <cmath> #include <cmath>
#include <cstring> #include <cstring>
#include <QDebug>
#include "Solve.hpp"
#include "Orbit.hpp"
using namespace std; using namespace std;
#define EPSILON 1e-10 #define EPSILON 1e-10
//#define EPSILON 1e-4
// line from vsop87.c, but also used by Heafner, 5.3.12. This is mu, we ign
ore the comet's mass i.r.t. the Sun's.
#define GAUSS_GRAV_CONST (0.01720209895*0.01720209895)
#if defined(_MSC_VER) #if defined(_MSC_VER)
// cuberoot is missing in VC++ !? // cuberoot is missing in VC++ !?
#define cbrt(x) pow((x),1./3.) #define cbrt(x) pow((x),1./3.)
#endif #endif
static //! Solve true anomaly nu for hyperbolic orbit.
void InitHyp(double q,double n,double e,double dt,double &a1,double &a2) { //! @param q: perihel distance
//! @param n: mean motion
//! @param e: excentricity
//! @param dt: days from perihel
//! @param rCosNu: r*cos(nu)
//! @param rSinNu: r*sin(nu)
static void InitHyp(const double q, const double n, const double e, const d
ouble dt, double &rCosNu, double &rSinNu) {
// qDebug() << "InitHyp";
Q_ASSERT(e>1.0);
const double a = q/(e-1.0); const double a = q/(e-1.0);
Q_ASSERT(a>0.0);
const double M = n * dt; const double M = n * dt;
double H = M; // double H = M;
for (;;) { // Newton // for (;;) { // Newton
const double Hp = H; // const double Hp = H;
H = H-(e*sinh(H)-H-M)/(e*cosh(H)-1); // H = H-(e*sinh(H)-H-M)/(e*cosh(H)-1);
if (fabs(H - Hp) < EPSILON) break; // if (fabs(H - Hp) < EPSILON) break;
// }
// const double h1 = q*sqrt((e+1.0)/(e-1.0));
// a1 = a*(e-cosh(H));
// a2 = h1*sinh(H);
// GZ Again I prefer Heafner, ch.5.4
double E=StelUtils::sign(M)*log(2.0*fabs(M)/e + 1.85);
// qDebug() << "InitHyp: E=" << E << " M=" << M ;
for (;;){
const double Ep=E;
const double f2=e*sinh(E);
const double f=f2-E-M;
const double f1=e*cosh(E)-1.0;
E+= (-5.0*f)/(f1+StelUtils::sign(f1)*sqrt(fabs(16.0*f1*f1-20.0*f*f2))
);
if (fabs(E-Ep) < EPSILON) break;
} }
const double h1 = q*sqrt((e+1.0)/(e-1.0)); rCosNu = a*(e-cosh(E));
a1 = a*(e-cosh(H)); rSinNu = a*sqrt(e*e-1.0)*sinh(E);
a2 = h1*sinh(H);
} }
static //! Solve true anomaly nu for parabolic orbit.
void InitPar(double q,double n,double dt,double &a1,double &a2) { //! @param q: perihel distance
//! @param n: mean motion equivalent related to W (n=W/dt) in Heafner, ch5.
5
//! @param dt: days from perihel
//! @param rCosNu: r*cos(nu)
//! @param rSinNu: r*sin(nu)
/*
static void InitPar(const double q, const double n, const double dt, doub
le &a1, double &a2) {
const double A = n*dt; const double A = n*dt;
const double h = sqrt(A*A+1.0); const double h = sqrt(A*A+1.0);
double c = cbrt(fabs(A)+h); double c = cbrt(fabs(A)+h);
c = c*c; c = c*c;
const double tan_nu_h = 2*A/(1+c+1/c); const double tan_nu_h = 2*A/(1+c+1/c);
a1 = q*(1-tan_nu_h*tan_nu_h); a1 = q*(1-tan_nu_h*tan_nu_h);
a2 = 2.0*q*tan_nu_h; a2 = 2.0*q*tan_nu_h;
} }
*/
static // GZ This implementation now follows Heafner, ch 5.5
void InitEll(double q,double n,double e,double dt,double &a1,double &a2) { static void InitPar(const double q, const double n, const double dt, double
const double a = q/(1.0-e); &rCosNu, double &rSinNu) {
double M = fmod(n*dt,2*M_PI); // qDebug() << "InitPar";
// const double M=dt*sqrt(GAUSS_GRAV_CONST/(2.0*q*q*q));
// const double W=1.5*M;
const double W=dt*n;
const double Y=cbrt(W+sqrt(W*W+1));
const double tanNu2=Y-1.0/Y; // Heafner (5.5.8) has an error here, writ
es (Y-1)/Y.
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)
//! @param q: perihel distance
//! @param n: mean motion
//! @param e: excentricity
//! @param dt: days from perihel
//! @param rCosNu: r*cos(nu)
//! @param rSinNu: r*sin(nu)
static void InitEll(const double q, const double n, const double e, const d
ouble dt, double &rCosNu, double &rSinNu) {
// qDebug() << "InitEll";
Q_ASSERT(e<1.0);
const double a = q/(1.0-e); // semimajor axis
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 H = M; // double E = M;
for (;;) { // Newton // for (;;) { // Newton(?) Solve Kepler's equation (similar to Meeus secon
const double Hp = H; d method, Astro.Alg. 1998 p.199)
H = H-(M-H+e*sin(H))/(e*cos(H)-1); // const double Ep = E;
if (fabs(H-Hp) < 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.
// Better always use Laguerre-Conway. See Heafner, Ch. 5.3
double E=M+0.85*e*StelUtils::sign(sin(M));
for (;;){
const double Ep=E;
const double f2=e*sin(E);
const double f=E-f2-M;
const double f1=1.0-e*cos(E);
E+= (-5.0*f)/(f1+StelUtils::sign(f1)*sqrt(fabs(16.0*f1*f1-20.0*f*f2))
);
if (fabs(E-Ep) < EPSILON) break;
} }
const double h1 = q*sqrt((1.0+e)/(1.0-e)); // Note: q=a*(1-e)
a1 = a*(cos(H)-e); const double h1 = q*sqrt((1.0+e)/(1.0-e)); // elsewhere: a sqrt(1-e²)
a2 = h1*sin(H); ... q / (1-e) sqrt( (1+e)(1-e)) = q sqrt((1+e)/(1-e))
} rCosNu = a*(cos(E)-e);
rSinNu = h1*sin(E);
void Init3D(double i,double Omega,double o,double a1,double a2, }
double &x1,double &x2,double &x3) {
const double co = cos(o); //! Compute position vector and (optional) speed vector from orbital elemen
const double so = sin(o); ts and true anomaly components. See e.g. Heafner, Fund.Eph.Comp.1999
//! @param i inclination
//! @param Omega, longitude of ascending node
//! @param w omega, argument of pericenter
//! @param rCosNu: r*cos(nu)
//! @param rSinNu: r*sin(nu)
//! @param rx: x component of position vector, AU
//! @param ry: y 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 e excentricity (required if withVelVector=true)
//! @param q perihel distance, AU (required if withVelVector=true)
//! @param rdotx: x component of velocity vector, AU/d
//! @param rdoty: y 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,
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 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 d11=-so*sOm*ci+co*cOm; const double Px=-sw*sOm*ci+cw*cOm; // Heafner, 5.3.1 Px
const double d12=-co*sOm*ci-so*cOm; const double Qx=-cw*sOm*ci-sw*cOm; // Heafner, 5.3.4 Qx
const double d21= so*cOm*ci+co*sOm; const double Py= sw*cOm*ci+cw*sOm; // Heafner, 5.3.2 Py
const double d22= co*cOm*ci-so*sOm; const double Qy= cw*cOm*ci-sw*sOm; // Heafner, 5.3.5 Qy
const double d31= so*si; const double Pz= sw*si; // Heafner, 5.3.3 Pz
const double d32= co*si; const double Qz= cw*si; // Heafner, 5.3.6 Qz
x1 = d11*a1+d12*a2; rx = Px*rCosNu+Qx*rSinNu; // Heafner, 5.3.18 r
x2 = d21*a1+d22*a2; ry = Py*rCosNu+Qy*rSinNu;
x3 = d31*a1+d32*a2; rz = Pz*rCosNu+Qz*rSinNu;
if (withVelVector) {
const double r=sqrt(rSinNu*rSinNu+rCosNu*rCosNu);
const double sinNu=rSinNu/r;
const double cosNu=rCosNu/r;
const double p=q*(1.0+e);
const double sqrtMuP=sqrt(GAUSS_GRAV_CONST/p);
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 meanMotion, double orbitGoodDays,
double parentRotObliquity, double meanMotion, // GZ: for parabolic
double parentRotAscendingnode, s, this is W/dt in Heafner's lettering
double parentRotJ2000Longitude) double parentRotObliquity, // GZ: I don't see a
:q(pericenterDistance),e(eccentricity),i(inclination), ny use for this, should be 0 for comets.
Om(ascendingNode),o(argOfPerhelion),t0(timeAtPerihelion), double parentRotAscendingnode, // GZ: I don't see a
n(meanMotion) { ny use for this, should be 0 for comets.
const double c_obl = cos(parentRotObliquity); double parentRotJ2000Longitude) // GZ: I don't see a
const double s_obl = sin(parentRotObliquity); ny use for this, should be 0 for comets.
const double c_nod = cos(parentRotAscendingnode); //)
const double s_nod = sin(parentRotAscendingnode); :q(pericenterDistance),e(eccentricity),i(inclination),
const double cj = cos(parentRotJ2000Longitude); Om(ascendingNode),w(argOfPerhelion),t0(timeAtPerihelion),
const double sj = sin(parentRotJ2000Longitude); n(meanMotion), updateTails(true), orbitGood(orbitGoodDays) {
// qDebug() << "CometOrbit::()";
rdot.set(0.0, 0.0, 0.0);
const double c_obl = cos(parentRotObliquity); // 1?
const double s_obl = sin(parentRotObliquity); // 0?
const double c_nod = cos(parentRotAscendingnode); // 1?
const double s_nod = sin(parentRotAscendingnode); // 0?
const double cj = cos(parentRotJ2000Longitude); // 1?
const double sj = sin(parentRotJ2000Longitude); // 0?
// 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; rotateToVsop87[0] = c_nod*cj-s_nod*c_obl*sj; // 1 in case of comet orbit
rotateToVsop87[1] = -c_nod*sj-s_nod*c_obl*cj; ing the sun, these names are misleading, however they do the right thing.
rotateToVsop87[2] = s_nod*s_obl; rotateToVsop87[1] = -c_nod*sj-s_nod*c_obl*cj; // 0 OK, this is NOT an ide
rotateToVsop87[3] = s_nod*cj+c_nod*c_obl*sj; ntity matrix...
rotateToVsop87[4] = -s_nod*sj+c_nod*c_obl*cj; rotateToVsop87[2] = s_nod*s_obl; // 0
rotateToVsop87[5] = -c_nod*s_obl; rotateToVsop87[3] = s_nod*cj+c_nod*c_obl*sj; // 0
rotateToVsop87[6] = s_obl*sj; rotateToVsop87[4] = -s_nod*sj+c_nod*c_obl*cj; // 1
rotateToVsop87[7] = s_obl*cj; rotateToVsop87[5] = -c_nod*s_obl; // 0
rotateToVsop87[8] = c_obl; rotateToVsop87[6] = s_obl*sj; // 0
rotateToVsop87[7] = s_obl*cj; // 0
rotateToVsop87[8] = c_obl; // 1
// qDebug() << "CometOrbit::()...done";
} }
void CometOrbit::positionAtTimevInVSOP87Coordinates(double JD,double *v) co nst { void CometOrbit::positionAtTimevInVSOP87Coordinates(double JD, double *v, b ool updateVelocityVector) {
JD -= t0; JD -= t0;
double a1,a2; 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);
else if (e > 1.0) InitHyp(q,n,e,JD,a1,a2); if (e < 1.0) InitEll(q,n,e,JD,rCosNu,rSinNu); // GZ: After solving with L
else InitPar(q,n,JD,a1,a2); aguerre-Conway, I dare to go for 1.0.
double p0,p1,p2; else if (e > 1.0) {
Init3D(i,Om,o,a1,a2,p0,p1,p2); // qDebug() << "Hyperbolic orbit for ecc=" << e << ", i=" << i << ", w
=" << w << ", Mean Motion n=" << n;
InitHyp(q,n,e,JD,rCosNu,rSinNu);
}
else InitPar(q,n,JD,rCosNu,rSinNu);
double p0,p1,p2, s0, s1, s2;
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) {
rdot.set(s0, s1, s2); // TODO: ROTATE?
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,
skipping to change at line 230 skipping to change at line 351
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
double operator()(double x) const // GZ: note&add Heafner's initial guess for E!
{ double operator()(double E) const
double s = ecc * sin(x); {
double c = ecc * cos(x); double s = ecc * sin(E);
double f = x - s - M; double c = ecc * cos(E);
double f = E - s - M;
double f1 = 1 - c; double f1 = 1 - c;
double f2 = s; double f2 = s;
x += -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 x; 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
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) ));
return x; return x;
} }
}; };
typedef pair<double, double> Solution; typedef pair<double, double> Solution;
double EllipticalOrbit::eccentricAnomaly(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, M), M, 5); Solution sol = solveIteration_fixed(SolveKeplerFunc1(eccentricity, M), M, 5);
skipping to change at line 300 skipping to change at line 422
{ {
// 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 calcuation.
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)
{ {
// Nearly 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.
// Anyhow, Comets use CometOrbit class.
return M; return M;
} }
else else
{ {
// Laguerre-Conway method for hyperbolic (ecc > 1) orbits. // Laguerre-Conway method for hyperbolic (ecc > 1) orbits.
double E = log(2 * M / eccentricity + 1.85); double E = log(2 * M / eccentricity + 1.85);
Solution sol = solveIteration_fixed(SolveKeplerLaguerreConwayHyp(ec centricity, M), E, 30); Solution sol = solveIteration_fixed(SolveKeplerLaguerreConwayHyp(ec centricity, M), E, 30);
return sol.first; return sol.first;
} }
} }
Vec3d EllipticalOrbit::positionAtE(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 * sqrt(1 - eccentricity * eccentricity) * -sin(E);
} }
else if (eccentricity > 1.0) else if (eccentricity > 1.0) // N.B. This is odd at least: elliptical m ust have ecc<1!
{ {
double a = pericenterDistance / (1.0 - eccentricity); double a = pericenterDistance / (1.0 - eccentricity);
x = -a * (eccentricity - cosh(E)); x = -a * (eccentricity - cosh(E));
z = -a * sqrt(eccentricity * eccentricity - 1) * -sinh(E); z = -a * sqrt(eccentricity * eccentricity - 1) * -sinh(E);
} }
else else
{ {
// TODO: Handle parabolic orbits // TODO: Handle parabolic orbits
x = 0.0; x = 0.0;
z = 0.0; z = 0.0;
} }
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(double t) const Vec3d EllipticalOrbit::positionAtTime(const double JD) const
{ {
t = t - epoch;
double meanMotion = 2.0 * M_PI / period; double meanMotion = 2.0 * M_PI / period;
double meanAnomaly = meanAnomalyAtEpoch + t * 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];
skipping to change at line 370 skipping to change at line 493
//} //}
//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(double JD, double* 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 87[2]*pos[2]; v[0] = rotateToVsop87[0]*pos[0] + rotateToVsop87[1]*pos[1] + rotateToVsop 87[2]*pos[2];
v[1] = rotateToVsop87[3]*pos[0] + rotateToVsop87[4]*pos[1] + rotateToVsop 87[5]*pos[2]; v[1] = rotateToVsop87[3]*pos[0] + rotateToVsop87[4]*pos[1] + rotateToVsop 87[5]*pos[2];
v[2] = rotateToVsop87[6]*pos[0] + rotateToVsop87[7]*pos[1] + rotateToVsop 87[8]*pos[2]; v[2] = rotateToVsop87[6]*pos[0] + rotateToVsop87[7]*pos[1] + rotateToVsop 87[8]*pos[2];
} }
double EllipticalOrbit::getPeriod() const double EllipticalOrbit::getPeriod() const
{ {
return period; return period;
 End of changes. 31 change blocks. 
91 lines changed or deleted 244 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/