#include <orsa_orbit.h>
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 |
Definition at line 101 of file orsa_orbit.h.
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.
Definition at line 337 of file orsa_orbit.cc.
References Orbit::Compute(), orsa::GetG(), Body::mass(), Orbit::mu, Body::position(), and Body::velocity().
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 }
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 }
void Compute | ( | const Body & | b, | |
const Body & | ref_b, | |||
const UniverseTypeAwareTime & | epoch_in | |||
) | [inline] |
Definition at line 125 of file orsa_orbit.h.
References OrbitWithEpoch::Compute(), and OrbitWithEpoch::epoch.
00125 { 00126 epoch = epoch_in; 00127 Orbit::Compute(b, ref_b); 00128 }
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 OrbitWithEpoch::epoch.
Referenced by OrbitWithEpoch::Compute(), 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 }
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(), orsa::sincos(), 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 }
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().
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 }
Reimplemented from Orbit.
Definition at line 131 of file orsa_orbit.h.
Referenced by OptimizedOrbitPositions::PropagatedOrbit(), OptimizedOrbitPositions::PropagatedPosVel(), 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 }
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(), orsa::least_sq_df(), orsa::least_sq_diff_f(), orsa::least_sq_f(), orsa::OrbitDifferentialCorrectionsLeastSquares(), 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(), Orbit::GetE(), orsa::least_sq_df(), orsa::least_sq_diff_f(), orsa::least_sq_f(), orsa::OrbitDifferentialCorrectionsLeastSquares(), 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().
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(), orsa::OrbitDifferentialCorrectionsLeastSquares(), 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(), orsa::least_sq_df(), orsa::least_sq_diff_f(), orsa::least_sq_f(), orsa::OrbitDifferentialCorrectionsLeastSquares(), 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(), Orbit::GetE(), orsa::least_sq_df(), orsa::least_sq_diff_f(), orsa::least_sq_f(), orsa::OrbitDifferentialCorrectionsLeastSquares(), 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(), orsa::least_sq_df(), orsa::least_sq_diff_f(), orsa::least_sq_f(), orsa::OrbitDifferentialCorrectionsLeastSquares(), 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(), orsa::least_sq_df(), orsa::least_sq_diff_f(), orsa::least_sq_f(), orsa::OrbitDifferentialCorrectionsLeastSquares(), 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().