OrbitWithEpoch Class Reference

#include <orsa_orbit.h>

Inheritance diagram for OrbitWithEpoch:

Inheritance graph
[legend]
Collaboration diagram for OrbitWithEpoch:

Collaboration graph
[legend]
List of all members.

Public Member Functions

 OrbitWithEpoch ()
 OrbitWithEpoch (const Orbit &orbit)
 OrbitWithEpoch (const OrbitWithEpoch &orbit)
void Compute (const Vector &dr, const Vector &dv, const double mu, const UniverseTypeAwareTime &epoch_in)
void Compute (const Body &b, const Body &ref_b, const UniverseTypeAwareTime &epoch_in)
void RelativePosVel (Vector &relative_position, Vector &relative_velocity) const
void RelativePosVel (Vector &relative_position, Vector &relative_velocity, const UniverseTypeAwareTime &epoch_in) const
double Period () const
double GetE () const
void Compute (const Vector &dr, const Vector &dv, const double mu)
void Compute (const Body &b, const Body &ref_b)

Public Attributes

UniverseTypeAwareTime epoch
double libration_angle
double a
double e
double i
double omega_node
double omega_pericenter
double M
double mu

Constructor & Destructor Documentation

OrbitWithEpoch  )  [inline]
 

Definition at line 103 of file orsa_orbit.h.

00103 : Orbit() { }

OrbitWithEpoch const Orbit orbit  )  [inline]
 

Definition at line 104 of file orsa_orbit.h.

00104 : Orbit(orbit) { }

OrbitWithEpoch const OrbitWithEpoch orbit  )  [inline]
 

Definition at line 105 of file orsa_orbit.h.

References OrbitWithEpoch::epoch.

00105 : Orbit(orbit), epoch(orbit.epoch) { }


Member Function Documentation

void Compute const Body b,
const Body ref_b
[inherited]
 

Definition at line 337 of file orsa_orbit.cc.

References Orbit::Compute(), orsa::GetG(), and Orbit::mu.

00337                                                       {
00338     
00339     const Vector dr = b.position() - ref_b.position();
00340     const Vector dv = b.velocity() - ref_b.velocity();
00341     
00342     const double mu = GetG()*(b.mass()+ref_b.mass());
00343     
00344     Compute(dr,dv,mu);
00345   }

Here is the call graph for this function:

void Compute const Vector dr,
const Vector dv,
const double  mu
[inherited]
 

Definition at line 347 of file orsa_orbit.cc.

References Orbit::a, orsa::cos(), Orbit::e, orsa::ExternalProduct(), Orbit::i, Vector::Length(), Vector::LengthSquared(), Orbit::M, Orbit::mu, Orbit::omega_node, Orbit::omega_pericenter, orsa::pi, orsa::secure_acos(), orsa::secure_atan2(), orsa::secure_log(), orsa::secure_sqrt(), orsa::sin(), orsa::tan(), orsa::twopi, Vector::x, Vector::y, and Vector::z.

Referenced by Orbit::Compute().

00347                                                                                                           {
00348     
00349     /////////////////////////////////////////////////////
00350     // This alghoritm is taken from the swift package, 
00351     // ORBEL_XV2EL.F file written by M. Duncan.
00352     /////////////////////////////////////////////////////
00353     
00354     mu = mu_in;
00355     
00356     const double tiny = 1.0e-100; // about 4.0e-15
00357     
00358     // internals
00359     double  face,cape,capf,tmpf,cw,sw,w,u;
00360     int ialpha = 0;
00361     
00362     // angular momentum
00363     Vector h = ExternalProduct(relative_position,relative_velocity);
00364     
00365     double h2 = h.LengthSquared();
00366     double hh = h.Length();
00367     
00368     // inclination
00369     i = secure_acos(h.z/hh);
00370     
00371     // Compute longitude of ascending node omega_node and the argument of latitude u
00372     // double fac = secure_sqrt(secure_pow(h.x,2)+secure_pow(h.y,2))/h2;
00373     double fac = secure_sqrt(h.x*h.x+h.y*h.y)/h2;
00374     
00375     if (fac < tiny) {
00376       omega_node = 0.0;
00377       u = secure_atan2(relative_position.y, relative_position.x);
00378       if ( fabs(i-pi) < 10.0 * tiny ) u = -u;
00379     } else {
00380       omega_node = secure_atan2(h.x,-h.y);
00381       u = secure_atan2(relative_position.z/sin(i), relative_position.x*cos(omega_node)+relative_position.y*sin(omega_node));
00382     }
00383     
00384     if (omega_node < 0.0) omega_node += twopi;
00385     if (u < 0.0) u += twopi;
00386     
00387     //  Compute the radius r and velocity squared v2, and the dot
00388     //  product rdotv, the energy per unit mass energy 
00389     double r  = relative_position.Length();
00390     double v2 = relative_velocity.LengthSquared();
00391     
00392     double vdotr  = relative_position*relative_velocity;
00393     
00394     double energy = v2/2.0 - mu/r;
00395     
00396     // Determine type of conic section and label it via ialpha
00397     if (fabs(energy*r/mu) < tiny) {
00398       ialpha = 0;
00399     } else {
00400       if (energy < 0)  ialpha = -1;
00401       if (energy > 0)  ialpha = +1;
00402     }
00403     
00404     // Depending on the conic type, determine the remaining elements 
00405     
00406     // ellipse 
00407     if (ialpha == -1) {
00408       
00409       a   = -mu/(2.0*energy); 
00410       
00411       fac = 1.0 - h2/(mu*a); 
00412       
00413       if (fac > tiny) {
00414         e = secure_sqrt(fac);
00415         face = (a-r)/(a*e);
00416         
00417         if (face > 1.0) {
00418           cape = 0.0;
00419         } else {
00420           if (face > -1.0)
00421             cape = secure_acos(face);
00422           else
00423             cape = pi;
00424         }
00425         
00426         if (vdotr < 0.0) cape = twopi - cape;
00427         cw = (cos(cape)-e)/(1.0-e*cos(cape));                  
00428         sw = secure_sqrt(1.0-e*e)*sin(cape)/(1.0-e*cos(cape));  
00429         w = secure_atan2(sw,cw);
00430         if (w < 0.0) w += twopi;
00431       } else {
00432         e = 0.0;
00433         w = u;
00434         cape = u;
00435       }
00436       
00437       M = cape - e*sin(cape);
00438       omega_pericenter = u - w;
00439       if (omega_pericenter < 0) omega_pericenter += twopi;
00440       omega_pericenter = fmod(omega_pericenter,twopi);
00441       
00442     }
00443     
00444     // hyperbola
00445     if (ialpha == 1) {
00446       a   = mu/(2.0*energy); 
00447       fac = h2/(mu*a); 
00448       
00449       if (fac > tiny) {
00450         e = secure_sqrt(1.0+fac);
00451         tmpf = (a+r)/(a*e);
00452         if (tmpf < 1.0) tmpf = 1.0;
00453         
00454         capf = secure_log(tmpf+secure_sqrt(tmpf*tmpf-1.0));
00455         
00456         if (vdotr < 0.0) capf = -capf;
00457         
00458         cw = (e-cosh(capf))/(e*cosh(capf)-1.0); 
00459         sw = secure_sqrt(e*e-1.0)*sinh(capf)/(e*cosh(capf)-1.0);
00460         w = secure_atan2(sw,cw);
00461         if (w < 0.0) w += twopi;
00462       } else {
00463         // we only get here if a hyperbola is essentially a parabola 
00464         // so we calculate e and w accordingly to avoid singularities
00465         e = 1.0;
00466         tmpf = h2/(2.0*mu); 
00467         w = secure_acos(2.0*tmpf/r - 1.0);
00468         if (vdotr < 0.0) w = twopi - w;
00469         tmpf = (a+r)/(a*e);
00470         capf = secure_log(tmpf+secure_sqrt(tmpf*tmpf-1.0));
00471       }
00472       
00473       M = e * sinh(capf) - capf; 
00474       omega_pericenter = u - w;
00475       if (omega_pericenter < 0) omega_pericenter += twopi;
00476       omega_pericenter = fmod(omega_pericenter,twopi);
00477     }
00478 
00479     // parabola
00480     //  NOTE - in this case we use "a" to mean pericentric distance
00481     if (ialpha == 0) {
00482       a = 0.5*h2/mu;
00483       e = 1.0;
00484       w = secure_acos(2.0*a/r -1.0);
00485       if (vdotr < 0.0) w = twopi - w;
00486       tmpf = tan(0.5*w);
00487       
00488       M = tmpf*(1.0+tmpf*tmpf/3.0);
00489       omega_pericenter = u - w;
00490       if (omega_pericenter < 0) omega_pericenter += twopi;
00491       omega_pericenter = fmod(omega_pericenter,twopi);
00492     }
00493     
00494   }

Here is the call graph for this function:

void Compute const Body b,
const Body ref_b,
const UniverseTypeAwareTime epoch_in
[inline]
 

Definition at line 125 of file orsa_orbit.h.

References orsa::Compute(), and OrbitWithEpoch::epoch.

00125                                                                                                     {
00126       epoch = epoch_in;
00127       Orbit::Compute(b, ref_b);
00128     }

Here is the call graph for this function:

void Compute const Vector dr,
const Vector dv,
const double  mu,
const UniverseTypeAwareTime epoch_in
[inline]
 

Definition at line 114 of file orsa_orbit.h.

References orsa::Compute(), and OrbitWithEpoch::epoch.

Referenced by orsa::Compute_Gauss(), orsa::Compute_TestMethod(), orsa::gauss_v_df(), orsa::gauss_v_diff_f(), orsa::gauss_v_f(), OptimizedOrbitPositions::PropagatedOrbit(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), and AstorbFile::Read().

00114                                                                                                                       {
00115       epoch = epoch_in;
00116       Orbit::Compute(dr,dv,mu);
00117     }

Here is the call graph for this function:

double GetE  )  const [inherited]
 

Definition at line 52 of file orsa_orbit.cc.

References orsa::cos(), orsa::E, Orbit::e, Orbit::M, ORSA_ERROR, ORSA_LOGIC_ERROR, ORSA_WARNING, orsa::pi, orsa::secure_pow(), orsa::sin(), and orsa::twopi.

Referenced by Orbit::RelativePosVel().

00052                            {
00053     
00054     if (e >= 1.0) {
00055       ORSA_WARNING("orsa::Orbit::GetE() called with eccentricity = %g; returning M.",e);
00056       //
00057       return M;
00058     }
00059 
00060     double E = 0.0;
00061     if (e < 0.8) {
00062       
00063 #ifdef HAVE_SINCOS
00064       double s,c;
00065       //
00066       ::sincos(M,&s,&c);
00067       const double sm = s;
00068       const double cm = c;
00069 #else // HAVE_SINCOS
00070       const double sm = sin(M);
00071       const double cm = cos(M);
00072 #endif // HAVE_SINCOS
00073       
00074       // begin with a guess accurate to order e^3 
00075       double x = M + e*sm*( 1.0 + e*( cm + e*( 1.0 -1.5*sm*sm)));
00076       
00077       double sx,cx;
00078       E = x;
00079       double old_E;
00080       double es,ec,f,fp,fpp,fppp,dx;
00081       
00082       unsigned int count = 0;
00083       const unsigned int max_count = 128;
00084       do {
00085         
00086 #ifdef HAVE_SINCOS
00087         ::sincos(x,&sx,&cx);
00088 #else // HAVE_SINCOS
00089         sx = sin(x);
00090         cx = cos(x);
00091 #endif // HAVE_SINCOS
00092         
00093         es = e*sx;
00094         ec = e*cx;
00095         f = x - es  - M;
00096         fp = 1.0 - ec; 
00097         fpp = es;
00098         fppp = ec; 
00099         dx = -f/fp;
00100         dx = -f/(fp + dx*fpp/2.0);
00101         dx = -f/(fp + dx*fpp/2.0 + dx*dx*fppp/6.0);
00102         //
00103         old_E = E;
00104         E = x + dx;
00105         ++count;
00106         // update x, ready for the next iteration
00107         x = E;
00108         
00109       } while ((fabs(E-old_E) > 10*(fabs(E)+fabs(M))*std::numeric_limits<double>::epsilon()) && (count < max_count));
00110       
00111       if (count >= max_count) {
00112         ORSA_ERROR("Orbit::GetE(): max count reached, e = %g    E = %g   fabs(E-old_E) = %g   10*(fabs(E)+fabs(M))*std::numeric_limits<double>::epsilon() = %g",e,E,fabs(E-old_E),10*(fabs(E)+fabs(M))*std::numeric_limits<double>::epsilon());
00113       }
00114            
00115     } else {
00116       
00117       double m = fmod(10*twopi+fmod(M,twopi),twopi);
00118       bool iflag = false;
00119       if (m > pi) {
00120         m = twopi - m;
00121         iflag = true;
00122       }
00123       
00124       // Make a first guess that works well for e near 1.  
00125       double x = secure_pow(6.0*m,1.0/3.0) - m;
00126       E = x;
00127       double old_E;
00128       double sa,ca,esa,eca,f,fp,dx;
00129        
00130       unsigned int count = 0;
00131       const unsigned int max_count = 128;
00132       do {
00133         
00134 #ifdef HAVE_SINCOS
00135         ::sincos(x+m,&sa,&ca);
00136 #else // HAVE_SINCOS
00137         sa = sin(x+m);
00138         ca = cos(x+m);
00139 #endif // HAVE_SINCOS
00140         
00141         esa = e*sa;
00142         eca = e*ca;
00143         f = x - esa;
00144         fp = 1.0 - eca;
00145         dx = -f/fp;
00146         dx = -f/(fp + 0.5*dx*esa);
00147         dx = -f/(fp + 0.5*dx*(esa+1.0/3.0*eca*dx));
00148         x += dx;
00149         //
00150         old_E = E;
00151         E = x + m;
00152         ++count;
00153         
00154       } while ((fabs(E-old_E) > 10*(fabs(E)+fabs(M)+twopi)*std::numeric_limits<double>::epsilon()) && (count < max_count));
00155       
00156       if (iflag) {
00157         E = twopi - E;
00158         old_E = twopi - old_E;
00159       }
00160       
00161       if (count >= max_count) {
00162         ORSA_WARNING("Orbit::GetE(): max count reached, e = %g    E = %g   fabs(E-old_E) = %g   10*(fabs(E)+fabs(M))*std::numeric_limits<double>::epsilon() = %g",e,E,fabs(E-old_E),10*(fabs(E)+fabs(M))*std::numeric_limits<double>::epsilon());
00163       }
00164     }
00165 
00166     if (E == 0.0) {
00167       ORSA_LOGIC_ERROR("E==0.0 in orsa::Orbit::GetE(); this should never happen.");
00168     }
00169     return E;
00170   }

Here is the call graph for this function:

double Period  )  const [inline, inherited]
 

Definition at line 84 of file orsa_orbit.h.

References Orbit::a, Orbit::mu, orsa::pisq, and orsa::secure_sqrt().

Referenced by OptimizedOrbitPositions::PropagatedOrbit(), JPLDastcomCometFile::Read(), MPCCometFile::Read(), and OrbitWithEpoch::RelativePosVel().

00084 { return (secure_sqrt(4*pisq*a*a*a/mu)); }

Here is the call graph for this function:

void RelativePosVel Vector relative_position,
Vector relative_velocity,
const UniverseTypeAwareTime epoch_in
const [inline]
 

Definition at line 136 of file orsa_orbit.h.

References OrbitWithEpoch::epoch, Orbit::M, Orbit::Period(), OrbitWithEpoch::RelativePosVel(), UniverseTypeAwareTime::Time(), and orsa::twopi.

00136                                                                                                                                     {
00137       OrbitWithEpoch o(*this);
00138       o.M += twopi*(epoch_in.Time() - epoch.Time())/Period();
00139       o.M  = std::fmod(10*twopi+std::fmod(o.M,twopi),twopi);
00140       o.RelativePosVel(relative_position,relative_velocity);
00141     }

Here is the call graph for this function:

void RelativePosVel Vector relative_position,
Vector relative_velocity
const [inline]
 

Reimplemented from Orbit.

Definition at line 131 of file orsa_orbit.h.

Referenced by OptimizedOrbitPositions::PropagatedOrbit(), OptimizedOrbitPositions::PropagatedPosVel(), OptimizedOrbitPositions::PropagatedSky_J2000(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), AstorbFile::Read(), and OrbitWithEpoch::RelativePosVel().

00131                                                                                              {
00132       Orbit::RelativePosVel(relative_position,relative_velocity);
00133     }


Member Data Documentation

double a [inherited]
 

Definition at line 97 of file orsa_orbit.h.

Referenced by Orbit::Compute(), orsa::Compute_Gauss(), orsa::Compute_TestMethod(), orsa::gauss_v_f(), OrbitWithCovarianceMatrixGSL::GenerateUsingCholeskyDecomposition(), OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation(), orsa::least_sq_df(), orsa::least_sq_f(), Orbit::Period(), Mercury5IntegrationFile::Read(), TLEFile::Read(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), AstorbFile::Read(), and Orbit::RelativePosVel().

double e [inherited]
 

Definition at line 97 of file orsa_orbit.h.

Referenced by Orbit::Compute(), orsa::Compute_Gauss(), orsa::Compute_TestMethod(), orsa::gauss_v_f(), OrbitWithCovarianceMatrixGSL::GenerateUsingCholeskyDecomposition(), OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation(), Orbit::GetE(), orsa::least_sq_df(), orsa::least_sq_f(), orsa::poly_8_gsl_solve(), Mercury5IntegrationFile::Read(), TLEFile::Read(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), AstorbFile::Read(), and Orbit::RelativePosVel().

UniverseTypeAwareTime epoch
 

Definition at line 144 of file orsa_orbit.h.

Referenced by OrbitWithEpoch::Compute(), orsa::gauss_v_df(), orsa::least_sq_df(), orsa::least_sq_f(), OrbitWithEpoch::OrbitWithEpoch(), OptimizedOrbitPositions::PropagatedOrbit(), Mercury5IntegrationFile::Read(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), AstorbFile::Read(), and OrbitWithEpoch::RelativePosVel().

double i [inherited]
 

Definition at line 97 of file orsa_orbit.h.

Referenced by Orbit::Compute(), orsa::Compute_Gauss(), orsa::Compute_TestMethod(), orsa::gauss_v_f(), OrbitWithCovarianceMatrixGSL::GenerateUsingCholeskyDecomposition(), OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation(), orsa::least_sq_df(), orsa::least_sq_f(), Mercury5IntegrationFile::Read(), TLEFile::Read(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), AstorbFile::Read(), and Orbit::RelativePosVel().

double libration_angle
 

Definition at line 145 of file orsa_orbit.h.

Referenced by SWIFTFile::Read().

double M [inherited]
 

Definition at line 97 of file orsa_orbit.h.

Referenced by Orbit::Compute(), OrbitWithCovarianceMatrixGSL::GenerateUsingCholeskyDecomposition(), OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation(), Orbit::GetE(), orsa::least_sq_df(), orsa::least_sq_f(), orsa::MOID(), orsa::MOID2RB(), OptimizedOrbitPositions::PropagatedOrbit(), Mercury5IntegrationFile::Read(), TLEFile::Read(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), AstorbFile::Read(), OrbitWithEpoch::RelativePosVel(), and Orbit::RelativePosVel().

double mu [inherited]
 

Definition at line 98 of file orsa_orbit.h.

Referenced by Orbit::Compute(), orsa::least_sq_df(), orsa::least_sq_f(), Orbit::Period(), TLEFile::Read(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), AstorbFile::Read(), and Orbit::RelativePosVel().

double omega_node [inherited]
 

Definition at line 97 of file orsa_orbit.h.

Referenced by Orbit::Compute(), OrbitWithCovarianceMatrixGSL::GenerateUsingCholeskyDecomposition(), OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation(), orsa::least_sq_df(), orsa::least_sq_f(), orsa::MOID(), orsa::MOID2RB(), Mercury5IntegrationFile::Read(), TLEFile::Read(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), AstorbFile::Read(), and Orbit::RelativePosVel().

double omega_pericenter [inherited]
 

Definition at line 97 of file orsa_orbit.h.

Referenced by Orbit::Compute(), OrbitWithCovarianceMatrixGSL::GenerateUsingCholeskyDecomposition(), OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation(), orsa::least_sq_df(), orsa::least_sq_f(), orsa::MOID(), orsa::MOID2RB(), Mercury5IntegrationFile::Read(), TLEFile::Read(), NEODYSCAT::Read(), JPLDastcomCometFile::Read(), JPLDastcomUnnumFile::Read(), JPLDastcomNumFile::Read(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), MPCCometFile::Read(), MPCOrbFile::Read(), AstorbFile::Read(), and Orbit::RelativePosVel().


The documentation for this class was generated from the following file:
Generated on Tue Jan 11 15:29:32 2005 for liborsa by  doxygen 1.4.0