orsa_frame.cc

Go to the documentation of this file.
00001 /* 
00002    ORSA - Orbit Reconstruction, Simulation and Analysis
00003    Copyright (C) 2002-2004 Pasquale Tricarico
00004    
00005    This program is free software; you can redistribute it and/or
00006    modify it under the terms of the GNU General Public License
00007    as published by the Free Software Foundation; either version 2
00008    of the License, or (at your option) any later version.
00009    
00010    As a special exception, Pasquale Tricarico gives permission to
00011    link this program with Qt commercial edition, and distribute the
00012    resulting executable, without including the source code for the Qt
00013    commercial edition in the source distribution.
00014    
00015    This program is distributed in the hope that it will be useful,
00016    but WITHOUT ANY WARRANTY; without even the implied warranty of
00017    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00018    GNU General Public License for more details.
00019    
00020    You should have received a copy of the GNU General Public License
00021    along with this program; if not, write to the Free Software
00022    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
00023 */
00024 
00025 #include "orsa_frame.h"
00026 
00027 #include "orsa_universe.h"
00028 #include "orsa_config.h"
00029 
00030 #include <iostream>
00031 
00032 using namespace std;
00033 
00034 namespace orsa {
00035   
00036   Frame::Frame() : vector<Body>(), UniverseTypeAwareTime() { }
00037   
00038   Frame::Frame(const UniverseTypeAwareTime & t) : vector<Body>(), UniverseTypeAwareTime(t) { }
00039   
00040   Frame::Frame(const Frame & f) : vector<Body>(f), UniverseTypeAwareTime(f) {
00041     
00042   }
00043   
00044   bool Frame::operator < (const Frame & f) const {
00045     return (this->UniverseTypeAwareTime::operator < (f));
00046   }
00047   
00048   void Frame::ForceJPLEphemerisData() {
00049     if (universe->GetUniverseType() != Real) return;
00050     for (unsigned int k=0;k<size();++k) {
00051       const JPL_planets p = (*this)[k].JPLPlanet();
00052       if (p == NONE) continue;
00053       JPLBody jb(p,*this);
00054       (*this)[k].SetPosition(jb.position());
00055       (*this)[k].SetVelocity(jb.velocity());      
00056     }
00057   }
00058   
00059   void print(const Frame &f) {
00060     cout << "Frame time: " << f.Time() << endl;
00061     cout << "Frame size: " << f.size() << endl;
00062     unsigned int l;
00063     for (l=0;l<f.size();l++) print(f[l]);
00064   }
00065   
00066   Vector Frame::CenterOfMass() const {
00067     return (orsa::CenterOfMass(*this));
00068   }
00069   
00070   Vector Frame::CenterOfMassVelocity() const {
00071     return (orsa::CenterOfMassVelocity(*this));
00072   }
00073   
00074   Vector Frame::Barycenter() const {
00075     return (orsa::Barycenter(*this));
00076   }
00077   
00078   Vector Frame::BarycenterVelocity() const {
00079     return (orsa::BarycenterVelocity(*this));
00080   }
00081   
00082   // Relativistic Barycenter stuff, from:
00083   // "Explanatory Supplement to the Astronomical Almanac",
00084   // Edited by P. Kenneth Seidelmann, U.S. Naval Observatory,
00085   // 2nd edition, 1992, sect. 5.215, pag. 286.
00086   
00087   // Eq. 5.215-2 of the book.
00088   // mu = g * mass
00089   static double modified_mu(const vector<Body> & f, const unsigned int i) {
00090     if (f[i].has_zero_mass()) return 0.0;
00091     const double one_over_two_c = 1.0/(2.0*GetC());
00092     // const double g = GetG();
00093     //
00094     vector<double> mu(f.size());
00095     for (unsigned int j=0; j<f.size(); ++j) {
00096       if (f[j].has_zero_mass()) {
00097         mu[j] = 0.0;
00098       } else {
00099         mu[j] = f[j].mu();
00100       }
00101     }
00102     //
00103     double mod_mu = 0.0;
00104     double tmp_sum = 0.0;
00105     for (unsigned int j=0; j<f.size(); ++j) {   
00106       if (j == i) continue;
00107       tmp_sum += mu[j]/(f[j].position()-f[i].position()).Length();
00108     }
00109     //
00110     mod_mu += mu[i]*(1.0+one_over_two_c*(f[i].velocity().LengthSquared()-tmp_sum));
00111     return mod_mu;
00112   }
00113   
00114   Vector Frame::RelativisticBarycenter() const {
00115     return (orsa::RelativisticBarycenter(*this));
00116   }
00117   
00118   Vector Frame::RelativisticBarycenterVelocity() const {
00119     return (orsa::RelativisticBarycenterVelocity(*this));
00120   }
00121   
00122   Vector CenterOfMass(const vector<Body> & f) {
00123     Vector sum_vec(0,0,0);
00124     double sum_mass = 0.0;
00125     for (unsigned int k=0; k<f.size(); ++k) {
00126       const double mass = f[k].mass();
00127       if (mass > 0.0) {
00128         sum_vec  += mass*f[k].position();
00129         sum_mass += mass;
00130       }
00131     }
00132     return (sum_vec/sum_mass);
00133   }
00134   
00135   Vector CenterOfMassVelocity(const vector<Body> & f) {
00136     Vector sum_vec(0,0,0);
00137     double sum_mass = 0.0;
00138     for (unsigned int k=0; k<f.size(); ++k) {
00139       const double mass = f[k].mass();
00140       if (mass > 0.0) {
00141         sum_vec  += mass*f[k].velocity();
00142         sum_mass += mass;
00143       }
00144     }
00145     return (sum_vec/sum_mass);
00146   }
00147   
00148   Vector Barycenter(const vector<Body> & f) {
00149     return (orsa::CenterOfMass(f));
00150   }
00151   
00152   Vector BarycenterVelocity(const vector<Body> & f) {
00153     return (orsa::CenterOfMassVelocity(f));
00154   }
00155   
00156   Vector RelativisticBarycenter(const vector<Body> & f) {
00157     Vector sum_vec(0,0,0);
00158     double sum_mu = 0.0;
00159     for (unsigned int i=0; i<f.size(); ++i) {
00160       const double mod_mu = modified_mu(f,i);
00161       if (mod_mu > 0.0) {
00162         sum_vec += mod_mu * f[i].position();
00163         sum_mu  += mod_mu;
00164       }
00165     }
00166     //
00167     return (sum_vec/sum_mu);
00168   }
00169   
00170   Vector RelativisticBarycenterVelocity(const vector<Body> & f) {
00171     Vector sum_vec(0,0,0);
00172     double sum_mu = 0.0;
00173     for (unsigned int i=0; i<f.size(); ++i) {
00174       const double mod_mu = modified_mu(f,i);
00175       if (mod_mu > 0.0) {
00176         sum_vec += mod_mu * f[i].velocity();
00177         sum_mu  += mod_mu;
00178       }
00179     }
00180     //
00181     return (sum_vec/sum_mu);
00182   }
00183   
00184   Vector Frame::AngularMomentum(const Vector & center) const {
00185     return (orsa::AngularMomentum(*this,center));
00186   }
00187   
00188   Vector AngularMomentum(const vector<Body> & f, const Vector & center) {
00189     Vector vec_sum;
00190     unsigned int k;
00191     for (k=0;k<f.size();++k) {
00192       if (f[k].mass() > 0) {
00193         vec_sum += f[k].mass()*ExternalProduct(f[k].position()-center,f[k].velocity());
00194       }
00195     }
00196     return (vec_sum);
00197   }
00198   
00199   Vector Frame::BarycentricAngularMomentum() const {
00200     return (orsa::BarycentricAngularMomentum(*this));
00201   }
00202   
00203   Vector BarycentricAngularMomentum(const vector<Body> & f) {
00204     return (orsa::AngularMomentum(f,Barycenter(f)));
00205   }
00206   
00207   double KineticEnergy(const vector<Body> & f) {
00208     if (f.size()==0) return (0.0);
00209     double energy=0.0;
00210     unsigned int k=0;
00211     while (k<f.size()) {
00212       energy += f[k].KineticEnergy();
00213       k++;
00214     }
00215     return (energy);
00216   }
00217   
00218 } // namespace orsa

Generated on Fri Nov 3 20:37:42 2006 for liborsa by  doxygen 1.4.7