PreliminaryOrbit Class Reference

#include <orsa_orbit_gsl.h>

Inheritance diagram for PreliminaryOrbit:

Inheritance graph
[legend]
Collaboration diagram for PreliminaryOrbit:

Collaboration graph
[legend]
List of all members.

Public Member Functions

double GetRMS () const
void ComputeRMS (const std::vector< Observation > &)
bool operator< (const PreliminaryOrbit &orbit) const
void SetCovarianceMatrix (double **, CovarianceMatrixElements)
void GetCovarianceMatrix (double **, CovarianceMatrixElements &) const
void GenerateUsingCholeskyDecomposition (std::vector< OrbitWithEpoch > &, const int, const int=12345) const
void GenerateUsingPrincipalAxisTransformation (std::vector< OrbitWithEpoch > &, const int, const int=12345) const
bool have_covariance_matrix () const
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 Compute (const Vector &dr, const Vector &dv, const double mu)
void Compute (const Body &b, const Body &ref_b)
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

Public Attributes

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

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, inherited]
 

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, inherited]
 

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:

void ComputeRMS const std::vector< Observation > &   ) 
 

Definition at line 63 of file orsa_orbit_gsl.cc.

References orsa::RMS_residuals().

00063                                                                   {
00064     rms = RMS_residuals(obs,*this);
00065   }

Here is the call graph for this function:

void GenerateUsingCholeskyDecomposition std::vector< OrbitWithEpoch > &  ,
const   int,
const   int = 12345
const [inherited]
 

Definition at line 2188 of file orsa_orbit_gsl.cc.

References Orbit::a, orsa::A, orsa::cos(), Orbit::e, orsa::Equinoctal, OrbitWithCovarianceMatrixGSL::have_covariance_matrix(), Orbit::i, Orbit::M, Orbit::omega_node, Orbit::omega_pericenter, ORSA_ERROR, ORSA_WARNING, orsa::Osculating, orsa::secure_atan2(), orsa::sin(), orsa::tan(), and orsa::twopi.

02188                                                                                                                                                  {
02189     
02190     list.clear();
02191     
02192     if (num < 1) return;
02193     
02194     if (!have_covariance_matrix()) {
02195       ORSA_ERROR("called method OrbitWithCovarianceMatrixGSL::GenerateUsingCholeskyDecomposition() from an orbit with undefined covariance matrix.");
02196       return;
02197     }
02198     
02199     gsl_vector * mu = gsl_vector_alloc(6);
02200     gsl_vector * x  = gsl_vector_alloc(6); // random vector
02201     gsl_vector * y  = gsl_vector_alloc(6);
02202     
02203     switch (cov_base) {
02204     case Osculating: {
02205       
02206       gsl_vector_set(mu,0,a);
02207       gsl_vector_set(mu,1,e);
02208       gsl_vector_set(mu,2,i);
02209       gsl_vector_set(mu,3,omega_node);
02210       gsl_vector_set(mu,4,omega_pericenter);
02211       gsl_vector_set(mu,5,M);
02212       
02213       break;
02214     }
02215     case Equinoctal: {
02216       
02217       gsl_vector_set(mu,0,a);
02218       gsl_vector_set(mu,1,e*sin(omega_node+omega_pericenter));
02219       gsl_vector_set(mu,2,e*cos(omega_node+omega_pericenter));
02220       gsl_vector_set(mu,3,tan(i/2)*sin(omega_node));
02221       gsl_vector_set(mu,4,tan(i/2)*cos(omega_node));
02222       gsl_vector_set(mu,5,fmod(omega_node+omega_pericenter+M,twopi));
02223       
02224       break;
02225     }
02226     }
02227     
02228     gsl_matrix * A = gsl_matrix_alloc(6,6);
02229     
02230     int i,k;
02231     for (i=0;i<6;++i) {
02232       for (k=0;k<6;++k) {
02233         gsl_matrix_set(A,i,k,covm[i][k]);
02234       }
02235     }
02236     
02237     // Cholesky decomposition
02238     int decomp_status = gsl_linalg_cholesky_decomp(A);
02239     
02240     if (decomp_status) {
02241       
02242       ORSA_WARNING("Cholesky decomposition failed.");
02243       
02244       gsl_vector_free(mu);
02245       gsl_vector_free(x);
02246       gsl_vector_free(y);
02247       gsl_matrix_free(A);
02248       
02249       return;
02250     }
02251     
02252     // remove L^T, the upper triangle (keep the diagonal elements!)
02253     for (i=0;i<6;++i) {
02254       for (k=i+1;k<6;++k) {
02255         gsl_matrix_set(A,i,k,0.0);
02256       }
02257     }
02258     
02259     // check
02260     // gsl_matrix_fprintf(stderr,A,"%g");
02261     
02262     OrbitWithEpoch gen_orb = (*this);
02263     
02264     int generated = 0;
02265     
02266     // random number generator
02267     gsl_rng * rnd = gsl_rng_alloc(gsl_rng_gfsr4);
02268     gsl_rng_set(rnd,random_seed);
02269     
02270     switch (cov_base) {
02271     case Osculating: {
02272       while (1) {
02273         
02274         for (i=0;i<6;++i) gsl_vector_set(x,i,gsl_ran_ugaussian(rnd));
02275         
02276         // y = A x + mu
02277         gsl_vector_memcpy(y,mu);
02278         gsl_blas_dgemv(CblasNoTrans,1.0,A,x,1.0,y);
02279         
02280         ++generated;
02281         
02282         gen_orb.a                 = gsl_vector_get (y,0);
02283         gen_orb.e                 = gsl_vector_get (y,1);
02284         gen_orb.i                 = gsl_vector_get (y,2);
02285         gen_orb.omega_node        = gsl_vector_get (y,3);
02286         gen_orb.omega_pericenter  = gsl_vector_get (y,4);
02287         gen_orb.M                 = gsl_vector_get (y,5);
02288         
02289         list.push_back(gen_orb);
02290         
02291         if (generated>=num) break;
02292       }
02293       break;
02294     }
02295     case Equinoctal: {
02296       
02297       while (1) {
02298         
02299         for (i=0;i<6;++i) gsl_vector_set(x,i,gsl_ran_ugaussian(rnd));
02300         
02301         // y = A x + mu
02302         gsl_vector_memcpy(y,mu);
02303         gsl_blas_dgemv(CblasNoTrans,1.0,A,x,1.0,y);
02304         
02305         ++generated;
02306         
02307         const double omega_tilde  = secure_atan2(gsl_vector_get(y,1),gsl_vector_get(y,2));
02308         
02309         gen_orb.a                 = gsl_vector_get(y,0);
02310         gen_orb.e                 = sqrt(gsl_vector_get(y,1)*gsl_vector_get(y,1)+gsl_vector_get(y,2)*gsl_vector_get(y,2));
02311         gen_orb.i                 = 2.0*atan(sqrt(gsl_vector_get(y,3)*gsl_vector_get(y,3)+gsl_vector_get(y,4)*gsl_vector_get(y,4)));
02312         gen_orb.omega_node        = fmod(10.0*twopi+secure_atan2(gsl_vector_get(y,3),gsl_vector_get(y,4)),twopi);
02313         gen_orb.omega_pericenter  = fmod(10.0*twopi+omega_tilde-omega_node,twopi);
02314         gen_orb.M                 = fmod(10.0*twopi+gsl_vector_get(y,5)-omega_tilde,twopi);
02315         
02316         list.push_back(gen_orb);
02317         
02318         if (generated>=num) break;
02319       }
02320       break;
02321     }
02322     }
02323     
02324     gsl_vector_free(mu);
02325     gsl_vector_free(x);
02326     gsl_vector_free(y);
02327     gsl_matrix_free(A);
02328     
02329     // free random number generator
02330     gsl_rng_free(rnd);
02331     
02332   }

Here is the call graph for this function:

void GenerateUsingPrincipalAxisTransformation std::vector< OrbitWithEpoch > &  ,
const   int,
const   int = 12345
const [inherited]
 

Definition at line 2336 of file orsa_orbit_gsl.cc.

References Orbit::a, orsa::A, orsa::cos(), Orbit::e, orsa::Equinoctal, OrbitWithCovarianceMatrixGSL::have_covariance_matrix(), Orbit::i, Orbit::M, Orbit::omega_node, Orbit::omega_pericenter, ORSA_ERROR, orsa::Osculating, orsa::secure_atan2(), orsa::secure_sqrt(), orsa::sin(), orsa::tan(), and orsa::twopi.

02336                                                                                                                                                        {
02337     
02338     list.clear();
02339     
02340     if (num < 1) return;
02341     
02342     if (!have_covariance_matrix()) {
02343       ORSA_ERROR("called method OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation() from an orbit with undefined covariance matrix.");
02344       return;
02345     }
02346     
02347     gsl_vector * mu = gsl_vector_alloc(6);
02348     gsl_vector * x  = gsl_vector_alloc(6); // random vector
02349     gsl_vector * y  = gsl_vector_alloc(6);
02350     
02351     switch (cov_base) {
02352     case Osculating: {
02353       
02354       gsl_vector_set(mu,0,a);
02355       gsl_vector_set(mu,1,e);
02356       gsl_vector_set(mu,2,i);
02357       gsl_vector_set(mu,3,omega_node);
02358       gsl_vector_set(mu,4,omega_pericenter);
02359       gsl_vector_set(mu,5,M);
02360       
02361       break;
02362     }
02363     case Equinoctal: {
02364       
02365       gsl_vector_set(mu,0,a);
02366       gsl_vector_set(mu,1,e*sin(omega_node+omega_pericenter));
02367       gsl_vector_set(mu,2,e*cos(omega_node+omega_pericenter));
02368       gsl_vector_set(mu,3,tan(i/2)*sin(omega_node));
02369       gsl_vector_set(mu,4,tan(i/2)*cos(omega_node));
02370       gsl_vector_set(mu,5,fmod(omega_node+omega_pericenter+M,twopi));
02371       
02372       break;
02373     }
02374     }
02375     
02376     gsl_matrix * A = gsl_matrix_alloc(6,6);
02377     
02378     for (unsigned int i=0;i<6;++i) {
02379       for (unsigned int k=0;k<6;++k) {
02380         gsl_matrix_set(A,i,k,covm[i][k]);
02381       }
02382     }
02383     
02384     // covariance matrix view
02385     // const gsl_matrix_view m = gsl_matrix_view_array(A,6,6);
02386     
02387     gsl_permutation * p = gsl_permutation_alloc(6);
02388     gsl_matrix * LU = gsl_matrix_alloc(6,6);
02389     gsl_matrix_memcpy(LU,A);  
02390     int s;
02391     
02392     // covariance matrix LU decomposition
02393     gsl_linalg_LU_decomp(LU, p, &s);
02394     
02395     // covariance matrix inverse
02396     gsl_matrix * inv = gsl_matrix_alloc(6,6);
02397     gsl_linalg_LU_invert(LU, p, inv);
02398     
02399     // covariance matrix determinant
02400     // const double det = gsl_linalg_LU_det(LU, s);
02401     
02402     // covariance eigenvalues and eigenvectors
02403     gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc(6); // workspace for eigenvectors/values
02404     gsl_vector * eval = gsl_vector_alloc(6);    // eigenvalues
02405     gsl_matrix * evec = gsl_matrix_alloc(6,6);  // eigenvectors
02406     //
02407     gsl_eigen_symmv (A, eval, evec, w); // NOTE: The diagonal and lower triangular part of A are destroyed during the computation
02408     
02409     /********/
02410     
02411     {
02412       // debug
02413       for (unsigned int i=0;i<6;++i) {
02414         cerr << "eval[" << i << "] = " << gsl_vector_get(eval,i) << endl;
02415       }
02416     }
02417     
02418     OrbitWithEpoch gen_orb = (*this);
02419     
02420     int generated = 0;
02421     
02422     // random number generator
02423     gsl_rng * rnd = gsl_rng_alloc(gsl_rng_gfsr4);
02424     gsl_rng_set(rnd,random_seed);
02425     
02426     double sigma[6];
02427     for (unsigned int i=0;i<6;++i) {
02428       if (gsl_vector_get(eval,i) < 0.0) {
02429         
02430         ORSA_ERROR("problems with the covariance matrix: negative eigenvalue found.");
02431         
02432         // free GSL stuff
02433         gsl_vector_free(mu);
02434         gsl_vector_free(x);
02435         gsl_vector_free(y);
02436         gsl_matrix_free(A);
02437         gsl_permutation_free(p);
02438         gsl_matrix_free(LU); 
02439         gsl_matrix_free(inv);
02440         gsl_eigen_symmv_free(w);
02441         gsl_vector_free(eval);
02442         gsl_matrix_free(evec);
02443         gsl_rng_free(rnd);
02444         
02445         return;
02446       }
02447       
02448       sigma[i] = secure_sqrt(gsl_vector_get(eval,i));
02449     }
02450     
02451     switch (cov_base) {
02452     case Osculating: {
02453       while (1) {
02454         
02455         for (unsigned int i=0;i<6;++i) {
02456           gsl_vector_set(x,i,gsl_ran_gaussian(rnd,sigma[i]));
02457         }
02458         
02459         gsl_vector_memcpy(y,mu);
02460         gsl_blas_dgemv(CblasNoTrans,1.0,evec,x,1.0,y);
02461         
02462         ++generated;
02463         
02464         gen_orb.a                 = gsl_vector_get (y,0);
02465         gen_orb.e                 = gsl_vector_get (y,1);
02466         gen_orb.i                 = gsl_vector_get (y,2);
02467         gen_orb.omega_node        = gsl_vector_get (y,3);
02468         gen_orb.omega_pericenter  = gsl_vector_get (y,4);
02469         gen_orb.M                 = gsl_vector_get (y,5);
02470         
02471         list.push_back(gen_orb);
02472         
02473         if (generated>=num) break;
02474       }
02475       break;
02476     }
02477     case Equinoctal: {
02478       while (1) {
02479         
02480         for (unsigned int i=0;i<6;++i) {
02481           gsl_vector_set(x,i,gsl_ran_gaussian(rnd,sigma[i]));
02482         }
02483         
02484         gsl_vector_memcpy(y,mu);
02485         gsl_blas_dgemv(CblasNoTrans,1.0,evec,x,1.0,y);
02486         
02487         ++generated;
02488         
02489         const double omega_tilde  = secure_atan2(gsl_vector_get(y,1),gsl_vector_get(y,2));
02490         
02491         gen_orb.a                 = gsl_vector_get(y,0);
02492         gen_orb.e                 = sqrt(gsl_vector_get(y,1)*gsl_vector_get(y,1)+gsl_vector_get(y,2)*gsl_vector_get(y,2));
02493         gen_orb.i                 = 2.0*atan(sqrt(gsl_vector_get(y,3)*gsl_vector_get(y,3)+gsl_vector_get(y,4)*gsl_vector_get(y,4)));
02494         gen_orb.omega_node        = fmod(10.0*twopi+secure_atan2(gsl_vector_get(y,3),gsl_vector_get(y,4)),twopi);
02495         gen_orb.omega_pericenter  = fmod(10.0*twopi+omega_tilde-omega_node,twopi);
02496         gen_orb.M                 = fmod(10.0*twopi+gsl_vector_get(y,5)-omega_tilde,twopi);
02497         
02498         list.push_back(gen_orb);
02499         
02500         if (generated>=num) break;
02501       }
02502       break;
02503     }
02504     }
02505     
02506     // free GSL stuff
02507     gsl_vector_free(mu);
02508     gsl_vector_free(x);
02509     gsl_vector_free(y);
02510     gsl_matrix_free(A);
02511     gsl_permutation_free(p);
02512     gsl_matrix_free(LU); 
02513     gsl_matrix_free(inv);
02514     gsl_eigen_symmv_free(w);
02515     gsl_vector_free(eval);
02516     gsl_matrix_free(evec);
02517     gsl_rng_free(rnd);
02518     
02519   } 

Here is the call graph for this function:

void GetCovarianceMatrix double **  ,
CovarianceMatrixElements
const [inherited]
 

Definition at line 2178 of file orsa_orbit_gsl.cc.

References OrbitWithCovarianceMatrixGSL::have_covariance_matrix(), and ORSA_ERROR.

02178                                                                                                                          {
02179     if (have_covariance_matrix()) {
02180       memcpy((void*)covariance_matrix, (const void*)covm, 6*6*sizeof(double));
02181       base = cov_base;
02182     } else {
02183       ORSA_ERROR("called method OrbitWithCovarianceMatrixGSL::GetCovarianceMatrix() from an orbit with undefined covariance matrix");
02184     }
02185   }

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 GetRMS  )  const [inline]
 

Definition at line 63 of file orsa_orbit_gsl.h.

00063 { return rms; };

bool have_covariance_matrix  )  const [inline, inherited]
 

Definition at line 55 of file orsa_orbit_gsl.h.

Referenced by OrbitWithCovarianceMatrixGSL::GenerateUsingCholeskyDecomposition(), OrbitWithCovarianceMatrixGSL::GenerateUsingPrincipalAxisTransformation(), and OrbitWithCovarianceMatrixGSL::GetCovarianceMatrix().

00055 { return bool_have_covariance_matrix; }

bool operator< const PreliminaryOrbit orbit  )  const [inline]
 

Definition at line 70 of file orsa_orbit_gsl.h.

00070                                                                  {
00071       return (rms < orbit.rms);
00072     }

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, inherited]
 

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, inherited]
 

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     }

void SetCovarianceMatrix double **  ,
CovarianceMatrixElements 
[inherited]
 

Definition at line 2172 of file orsa_orbit_gsl.cc.

Referenced by OrbitWithCovarianceMatrixGSL::OrbitWithCovarianceMatrixGSL().

02172                                                                                                                   {
02173     memcpy((void*)covm, (const void*)covariance_matrix, 6*6*sizeof(double));
02174     cov_base = base;
02175     bool_have_covariance_matrix = true;
02176   }


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 [inherited]
 

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 [inherited]
 

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 files:
Generated on Tue Jan 11 15:29:41 2005 for liborsa by  doxygen 1.4.0