00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
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
00083
00084
00085
00086
00087
00088
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
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 }