#include <orsa_orbit.h>
Inheritance diagram for Orbit:
Public Member Functions | |
Orbit () | |
double | Period () const |
double | GetE () const |
void | RelativePosVel (Vector &relative_position, Vector &relative_velocity) const |
void | Compute (const Vector &dr, const Vector &dv, const double mu) |
void | Compute (const Body &b, const Body &ref_b) |
Public Attributes | |
double | a |
double | e |
double | i |
double | omega_node |
double | omega_pericenter |
double | M |
double | mu |
Definition at line 79 of file orsa_orbit.h.
Orbit | ( | ) | [inline] |
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 }
Here is the call graph for this function:
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:
double GetE | ( | ) | const |
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 }
Here is the call graph for this function:
double Period | ( | ) | const [inline] |
Definition at line 84 of file orsa_orbit.h.
References Orbit::a, Orbit::mu, orsa::pisq, and orsa::secure_sqrt().
Referenced by OrbitWithEpoch::RelativePosVel().
Here is the call graph for this function:
Reimplemented in OrbitWithEpoch.
Definition at line 173 of file orsa_orbit.cc.
References Orbit::a, orsa::cos(), Orbit::e, Orbit::GetE(), Orbit::i, Orbit::M, Orbit::mu, Orbit::omega_node, Orbit::omega_pericenter, orsa::secure_log(), orsa::secure_sqrt(), orsa::sin(), orsa::sincos(), and orsa::twopi.
Referenced by orsa::MOID(), orsa::MOID2RB(), TLEFile::Read(), and OrbitWithEpoch::RelativePosVel().
00173 { 00174 00175 ///////////////////////////////////////////////////// 00176 // This alghoritm is taken from the swift package, 00177 // ORBEL_EL2XV.F file written by M. Duncan. 00178 ///////////////////////////////////////////////////// 00179 00180 // Generate rotation matrices (on p. 42 of Fitzpatrick) 00181 00182 // double sp,cp,so,co,si,ci; 00183 00184 #ifdef HAVE_SINCOS 00185 double s,c; 00186 00187 ::sincos(omega_pericenter,&s,&c); 00188 const double sp = s; 00189 const double cp = c; 00190 00191 ::sincos(omega_node,&s,&c); 00192 const double so = s; 00193 const double co = c; 00194 00195 ::sincos(i,&s,&c); 00196 const double si = s; 00197 const double ci = c; 00198 #else // HAVE_SINCOS 00199 const double sp = sin(omega_pericenter); 00200 const double cp = cos(omega_pericenter); 00201 00202 const double so = sin(omega_node); 00203 const double co = cos(omega_node); 00204 00205 const double si = sin(i); 00206 const double ci = cos(i); 00207 #endif // HAVE_SINCOS 00208 00209 const Vector d1(cp*co - sp*so*ci, 00210 cp*so + sp*co*ci, 00211 sp*si); 00212 00213 const Vector d2(-sp*co - cp*so*ci, 00214 -sp*so + cp*co*ci, 00215 cp*si); 00216 00217 // Get the other quantities depending on orbit type 00218 00219 // double cape,scap,ccap,sqe,sqgma,tmp; 00220 double cape,tmp; 00221 double xfac1,xfac2,vfac1,vfac2; 00222 00223 double capf,shcap,chcap; 00224 double zpara; 00225 00226 if (e < 1.0) { 00227 00228 /* 00229 int count = 0; 00230 cape = M; 00231 do { 00232 tmp = cape; 00233 cape = e * sin(cape) + M; 00234 ++count; 00235 } while ( (fabs((cape-tmp)/cape) > 1.0e-15) && (count < 100) ); 00236 */ 00237 00238 cape = GetE(); 00239 00240 // cerr << "tmp: " << tmp << " cape: " << cape << " fabs(cape - tmp)" << fabs(cape - tmp) << endl; 00241 00242 #ifdef HAVE_SINCOS 00243 ::sincos(cape,&s,&c); 00244 const double scap = s; 00245 const double ccap = c; 00246 #else // HAVE_SINCOS 00247 const double scap = sin(cape); 00248 const double ccap = cos(cape); 00249 #endif // HAVE_SINCOS 00250 const double sqe = secure_sqrt(1.0 - e*e); 00251 const double sqgma = secure_sqrt(mu*a); 00252 00253 xfac1 = a*(ccap - e); 00254 xfac2 = a*sqe*scap; 00255 // ri = 1/r 00256 const double ri = 1.0/(a*(1.0 - e*ccap)); 00257 vfac1 = -ri * sqgma * scap; 00258 vfac2 = ri * sqgma * ccap * sqe; 00259 00260 } else if (e > 1.0) { 00261 00262 double x,shx,chx,esh,ech; 00263 double f,fp,fpp,fppp,dx; 00264 00265 // use the 'right' value for M -- NEEDED! 00266 double local_M = M; 00267 if (fabs(local_M-twopi) < fabs(local_M)) { 00268 local_M -= twopi; 00269 } 00270 00271 // begin with a guess proposed by Danby 00272 if (local_M < 0.0) { 00273 tmp = -2.0*local_M/e + 1.8; 00274 x = -secure_log(tmp); 00275 } else { 00276 tmp = +2.0*local_M/e + 1.8; 00277 x = secure_log(tmp); 00278 } 00279 00280 capf = x; 00281 00282 int count = 0; 00283 do { 00284 x = capf; 00285 shx = sinh(x); 00286 chx = cosh(x); 00287 esh = e*shx; 00288 ech = e*chx; 00289 f = esh-x-local_M; 00290 fp = ech - 1.0; 00291 fpp = esh; 00292 fppp = ech; 00293 dx = -f/fp; 00294 dx = -f/(fp + dx*fpp/2.0); 00295 dx = -f/(fp + dx*fpp/2.0 + dx*dx*fppp/6.0); 00296 capf = x + dx; 00297 ++count; 00298 } while ((fabs(dx) > 1.0e-14) && (count < 100)); 00299 00300 shcap = sinh(capf); 00301 chcap = cosh(capf); 00302 00303 const double sqe = secure_sqrt(e*e - 1.0); 00304 const double sqgma = secure_sqrt(mu*a); 00305 xfac1 = a*(e - chcap); 00306 xfac2 = a*sqe*shcap; 00307 const double ri = 1.0/(a*(e*chcap - 1.0)); 00308 vfac1 = -ri * sqgma * shcap; 00309 vfac2 = ri * sqgma * chcap * sqe; 00310 00311 } else { // e = 1.0 within roundoff errors 00312 00313 double q = M; 00314 if (q < 1.0e-3) { 00315 zpara = q*(1.0 - (q*q/3.0)*(1.0-q*q)); 00316 } else { 00317 double x = 0.5*(3.0*q+secure_sqrt(9.0*(q*q)+4.0)); 00318 // double tmp = secure_pow(x,(1.0/3.0)); 00319 double tmp = cbrt(x); 00320 zpara = tmp - 1.0/tmp; 00321 } 00322 00323 const double sqgma = secure_sqrt(2.0*mu*a); 00324 xfac1 = a*(1.0 - zpara*zpara); 00325 xfac2 = 2.0*a*zpara; 00326 const double ri = 1.0/(a*(1.0 + zpara*zpara)); 00327 vfac1 = -ri * sqgma * zpara; 00328 vfac2 = ri * sqgma; 00329 00330 } 00331 00332 relative_position = d1*xfac1 + d2*xfac2; 00333 relative_velocity = d1*vfac1 + d2*vfac2; 00334 00335 }
Here is the call graph for this function:
double a |
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(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), and Orbit::RelativePosVel().
double e |
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(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), and Orbit::RelativePosVel().
double i |
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(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), and Orbit::RelativePosVel().
double M |
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::MOID(), orsa::MOID2RB(), orsa::OrbitDifferentialCorrectionsLeastSquares(), OptimizedOrbitPositions::PropagatedOrbit(), Mercury5IntegrationFile::Read(), TLEFile::Read(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), OrbitWithEpoch::RelativePosVel(), and Orbit::RelativePosVel().
double mu |
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(), AstDySMatrixFile::Read(), and Orbit::RelativePosVel().
double omega_node |
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::MOID(), orsa::MOID2RB(), orsa::OrbitDifferentialCorrectionsLeastSquares(), Mercury5IntegrationFile::Read(), TLEFile::Read(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), and Orbit::RelativePosVel().
double omega_pericenter |
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::MOID(), orsa::MOID2RB(), orsa::OrbitDifferentialCorrectionsLeastSquares(), Mercury5IntegrationFile::Read(), TLEFile::Read(), AstDySMatrixFile::Read(), RadauModIntegrationFile::Read(), SWIFTFile::Read(), and Orbit::RelativePosVel().