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_body.h"
00026 #include "orsa_error.h"
00027 #include "orsa_file.h"
00028
00029 #include <algorithm>
00030
00031 #include <iostream>
00032
00033 using namespace std;
00034
00035 namespace orsa {
00036
00037 unsigned int BodyConstants::used_body_id = 0;
00038
00039 list<BodyConstants*> BodyConstants::list_bc;
00040
00041 BodyConstants::BodyConstants() : name_(""), mass_(0.0), mu_(0.0), zero_mass_(mass_ == 0.0), radius_(0.0), planet_(NONE), J2_(0.0), J3_(0.0), J4_(0.0), C22_(0.0), C31_(0.0), C32_(0.0), C33_(0.0), C41_(0.0), C42_(0.0), C43_(0.0), C44_(0.0), S31_(0.0), S32_(0.0), S33_(0.0), S41_(0.0), S42_(0.0), S43_(0.0), S44_(0.0), id(used_body_id++) {
00042 users = 1;
00043 list_bc.push_back(this);
00044 }
00045
00046 BodyConstants::BodyConstants(const string & name, const double mass) : name_(name), mass_(mass), mu_(GetG()*mass_), zero_mass_(mass_ == 0.0), radius_(0.0), planet_(NONE), J2_(0.0), J3_(0.0), J4_(0.0), C22_(0.0), C31_(0.0), C32_(0.0), C33_(0.0), C41_(0.0), C42_(0.0), C43_(0.0), C44_(0.0), S31_(0.0), S32_(0.0), S33_(0.0), S41_(0.0), S42_(0.0), S43_(0.0), S44_(0.0), id(used_body_id++) {
00047 users = 1;
00048 list_bc.push_back(this);
00049 }
00050
00051 BodyConstants::BodyConstants(const string & name, const double mass, const double radius) : name_(name), mass_(mass), mu_(GetG()*mass_), zero_mass_(mass_ == 0.0), radius_(radius), planet_(NONE), J2_(0.0), J3_(0.0), J4_(0.0), C22_(0.0), C31_(0.0), C32_(0.0), C33_(0.0), C41_(0.0), C42_(0.0), C43_(0.0), C44_(0.0), S31_(0.0), S32_(0.0), S33_(0.0), S41_(0.0), S42_(0.0), S43_(0.0), S44_(0.0), id(used_body_id++) {
00052 users = 1;
00053 list_bc.push_back(this);
00054 }
00055
00056 BodyConstants::BodyConstants(const string & name, const double mass, const double radius, const JPL_planets p) : name_(name), mass_(mass), mu_(GetG()*mass_), zero_mass_(mass_ == 0.0), radius_(radius), planet_(p), J2_(0.0), J3_(0.0), J4_(0.0), C22_(0.0), C31_(0.0), C32_(0.0), C33_(0.0), C41_(0.0), C42_(0.0), C43_(0.0), C44_(0.0), S31_(0.0), S32_(0.0), S33_(0.0), S41_(0.0), S42_(0.0), S43_(0.0), S44_(0.0), id(used_body_id++) {
00057 users = 1;
00058 list_bc.push_back(this);
00059 }
00060
00061 BodyConstants::BodyConstants(const string & name, const double mass, const double radius, const double J2, const double J3, const double J4) : name_(name), mass_(mass), mu_(GetG()*mass_), zero_mass_(mass_ == 0.0), radius_(radius), planet_(NONE), J2_(J2), J3_(J3), J4_(J4), C22_(0.0), C31_(0.0), C32_(0.0), C33_(0.0), C41_(0.0), C42_(0.0), C43_(0.0), C44_(0.0), S31_(0.0), S32_(0.0), S33_(0.0), S41_(0.0), S42_(0.0), S43_(0.0), S44_(0.0), id(used_body_id++) {
00062 users = 1;
00063 list_bc.push_back(this);
00064 }
00065
00066 BodyConstants::BodyConstants(const string & name, const double mass, const double radius, const JPL_planets p, const double J2, const double J3, const double J4) : name_(name), mass_(mass), mu_(GetG()*mass_), zero_mass_(mass_ == 0.0), radius_(radius), planet_(p), J2_(J2), J3_(J3), J4_(J4), C22_(0.0), C31_(0.0), C32_(0.0), C33_(0.0), C41_(0.0), C42_(0.0), C43_(0.0), C44_(0.0), S31_(0.0), S32_(0.0), S33_(0.0), S41_(0.0), S42_(0.0), S43_(0.0), S44_(0.0), id(used_body_id++) {
00067 users = 1;
00068 list_bc.push_back(this);
00069 }
00070
00071 BodyConstants::BodyConstants(const string & name, const double mass, const double radius, const double J2, const double J3, const double J4, const double C22, const double C31, const double C32, const double C33, const double C41, const double C42, const double C43, const double C44, const double S31, const double S32, const double S33, const double S41, const double S42, const double S43, const double S44) : name_(name), mass_(mass), mu_(GetG()*mass_), zero_mass_(mass_ == 0.0), radius_(radius), planet_(NONE), J2_(J2), J3_(J3), J4_(J4), C22_(C22), C31_(C31), C32_(C32), C33_(C33), C41_(C41), C42_(C42), C43_(C43), C44_(C44), S31_(S31), S32_(S32), S33_(S33), S41_(S41), S42_(S42), S43_(S43), S44_(S44), id(used_body_id++) {
00072 users = 1;
00073 list_bc.push_back(this);
00074 }
00075
00076 BodyConstants::BodyConstants(const string & name, const double mass, const double radius, const JPL_planets p, const double J2, const double J3, const double J4, const double C22, const double C31, const double C32, const double C33, const double C41, const double C42, const double C43, const double C44, const double S31, const double S32, const double S33, const double S41, const double S42, const double S43, const double S44) : name_(name), mass_(mass), mu_(GetG()*mass_), zero_mass_(mass_ == 0.0), radius_(radius), planet_(p), J2_(J2), J3_(J3), J4_(J4), C22_(C22), C31_(C31), C32_(C32), C33_(C33), C41_(C41), C42_(C42), C43_(C43), C44_(C44), S31_(S31), S32_(S32), S33_(S33), S41_(S41), S42_(S42), S43_(S43), S44_(S44), id(used_body_id++) {
00077 users = 1;
00078 list_bc.push_back(this);
00079 }
00080
00081 BodyConstants::~BodyConstants() {
00082 list<BodyConstants*>::iterator it = list_bc.begin();
00083 while (it != list_bc.end()) {
00084 if (*it == this) {
00085 list_bc.erase(it);
00086 break;
00087 }
00088 ++it;
00089 }
00090 }
00091
00092
00093
00094 Body::Body() {
00095 bc = new BodyConstants();
00096 }
00097
00098 Body::Body(const double mass) {
00099 bc = new BodyConstants("",mass,0);
00100 }
00101
00102 Body::Body(const string & name) {
00103 bc = new BodyConstants(name,0,0);
00104 }
00105
00106 Body::Body(const string & name, const double mass) {
00107 bc = new BodyConstants(name,mass,0);
00108 }
00109
00110 Body::Body(const string & name, const double mass, const double radius) {
00111 bc = new BodyConstants(name,mass,radius);
00112 }
00113
00114 Body::Body(const string & name, const double mass, const double radius, const JPL_planets p) {
00115 bc = new BodyConstants(name,mass,radius,p);
00116 }
00117
00118 Body::Body(const string & name, const double mass, const Vector &position, const Vector &velocity) {
00119 bc = new BodyConstants(name,mass,0);
00120 _position = position;
00121 _velocity = velocity;
00122 }
00123
00124 Body::Body(const string & name, const double mass, const double radius, const Vector &position, const Vector &velocity) {
00125 bc = new BodyConstants(name,mass,radius);
00126 _position = position;
00127 _velocity = velocity;
00128 }
00129
00130 Body::Body(const string & name, const double mass, const double radius, const double J2, const double J3, const double J4) {
00131 bc = new BodyConstants(name,mass,radius,J2,J3,J4);
00132 }
00133
00134 Body::Body(const string & name, const double mass, const double radius, const JPL_planets p, const double J2, const double J3, const double J4) {
00135 bc = new BodyConstants(name,mass,radius,p,J2,J3,J4);
00136 }
00137
00138 Body::Body(const string & name, const double mass, const double radius, const double J2, const double J3, const double J4, const double C22, const double C31, const double C32, const double C33, const double C41, const double C42, const double C43, const double C44, const double S31, const double S32, const double S33, const double S41, const double S42, const double S43, const double S44) {
00139 bc = new BodyConstants(name,mass,radius,J2,J3,J4,C22,C31,C32,C33,C41,C42,C43,C44,S31,S32,S33,S41,S42,S43,S44);
00140 }
00141
00142 Body::Body(const string & name, const double mass, const double radius, const JPL_planets p, const double J2, const double J3, const double J4, const double C22, const double C31, const double C32, const double C33, const double C41, const double C42, const double C43, const double C44, const double S31, const double S32, const double S33, const double S41, const double S42, const double S43, const double S44) {
00143 bc = new BodyConstants(name,mass,radius,p,J2,J3,J4,C22,C31,C32,C33,C41,C42,C43,C44,S31,S32,S33,S41,S42,S43,S44);
00144 }
00145
00146 Body::Body(const string & name, const double mass, const Vector &position, const Vector &velocity, const double J2, const double J3, const double J4) {
00147 bc = new BodyConstants(name,mass,0,J2,J3,J4);
00148 _position = position;
00149 _velocity = velocity;
00150 }
00151
00152 Body::Body(const string & name, const double mass, const double radius, const Vector &position, const Vector &velocity, const double J2, const double J3, const double J4) {
00153 bc = new BodyConstants(name,mass,radius,J2,J3,J4);
00154 _position = position;
00155 _velocity = velocity;
00156 }
00157
00158 Body::Body(const Body & b) {
00159 bc = b.bc;
00160 bc->AddUser();
00161 _position = b._position;
00162 _velocity = b._velocity;
00163 }
00164
00165 Body::Body(const BodyWithEpoch & b) {
00166 bc = b.bc;
00167 bc->AddUser();
00168 _position = b._position;
00169 _velocity = b._velocity;
00170 }
00171
00172 Body::Body(const JPLBody & b) {
00173 bc = b.bc;
00174 bc->AddUser();
00175 _position = b.position();
00176 _velocity = b.velocity();
00177 }
00178
00179 Body::~Body() {
00180 bc->RemoveUser();
00181 if (bc->Users() == 0) {
00182 delete bc;
00183 bc = 0;
00184 }
00185 }
00186
00187 Body & Body::operator = (const Body & b) {
00188
00189 if (bc != b.bc) {
00190 bc->RemoveUser();
00191 if (bc->Users() == 0) {
00192 delete bc;
00193 bc = 0;
00194 }
00195
00196 bc = b.bc;
00197 bc->AddUser();
00198 }
00199
00200 _position = b._position;
00201 _velocity = b._velocity;
00202
00203 return * this;
00204 }
00205
00206
00207
00208 BodyWithEpoch::BodyWithEpoch(const BodyWithEpoch & b) : Body(b), epoch(b.epoch) { }
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243 static double local_mass(const JPL_planets planet) {
00244 return jpl_file->GetMass(planet);
00245 }
00246
00247 static double local_J2(const JPL_planets planet) {
00248 double J2_ = 0.0;
00249 switch (planet) {
00250 case SUN: J2_ = jpl_file->GetTag("J2SUN"); break;
00251 case EARTH: J2_ = jpl_file->GetTag("J2E"); break;
00252 case MOON: J2_ = jpl_file->GetTag("J2M"); break;
00253 default: break;
00254 }
00255 return J2_;
00256 }
00257
00258 static double local_J3(const JPL_planets planet) {
00259 double J3_ = 0.0;
00260 switch (planet) {
00261 case EARTH: J3_ = jpl_file->GetTag("J3E"); break;
00262 case MOON: J3_ = jpl_file->GetTag("J3M"); break;
00263 default: break;
00264 }
00265 return J3_;
00266 }
00267
00268 static double local_J4(const JPL_planets planet) {
00269 double J4_ = 0.0;
00270 switch (planet) {
00271 case EARTH: J4_ = jpl_file->GetTag("J4E"); break;
00272 case MOON: J4_ = jpl_file->GetTag("J4M"); break;
00273 default: break;
00274 }
00275 return J4_;
00276 }
00277
00278 static double local_C22(const JPL_planets planet) {
00279 if (planet == MOON) {
00280 return jpl_file->GetTag("C22M");
00281 } else {
00282 return 0.0;
00283 }
00284 }
00285
00286 static double local_C31(const JPL_planets planet) {
00287 if (planet == MOON) {
00288 return jpl_file->GetTag("C31M");
00289 } else {
00290 return 0.0;
00291 }
00292 }
00293
00294 static double local_C32(const JPL_planets planet) {
00295 if (planet == MOON) {
00296 return jpl_file->GetTag("C32M");
00297 } else {
00298 return 0.0;
00299 }
00300 }
00301
00302 static double local_C33(const JPL_planets planet) {
00303 if (planet == MOON) {
00304 return jpl_file->GetTag("C33M");
00305 } else {
00306 return 0.0;
00307 }
00308 }
00309
00310 static double local_C41(const JPL_planets planet) {
00311 if (planet == MOON) {
00312 return jpl_file->GetTag("C41M");
00313 } else {
00314 return 0.0;
00315 }
00316 }
00317
00318 static double local_C42(const JPL_planets planet) {
00319 if (planet == MOON) {
00320 return jpl_file->GetTag("C42M");
00321 } else {
00322 return 0.0;
00323 }
00324 }
00325
00326 static double local_C43(const JPL_planets planet) {
00327 if (planet == MOON) {
00328 return jpl_file->GetTag("C43M");
00329 } else {
00330 return 0.0;
00331 }
00332 }
00333
00334 static double local_C44(const JPL_planets planet) {
00335 if (planet == MOON) {
00336 return jpl_file->GetTag("C44M");
00337 } else {
00338 return 0.0;
00339 }
00340 }
00341
00342 static double local_S31(const JPL_planets planet) {
00343 if (planet == MOON) {
00344 return jpl_file->GetTag("S31M");
00345 } else {
00346 return 0.0;
00347 }
00348 }
00349
00350 static double local_S32(const JPL_planets planet) {
00351 if (planet == MOON) {
00352 return jpl_file->GetTag("S32M");
00353 } else {
00354 return 0.0;
00355 }
00356 }
00357
00358 static double local_S33(const JPL_planets planet) {
00359 if (planet == MOON) {
00360 return jpl_file->GetTag("S33M");
00361 } else {
00362 return 0.0;
00363 }
00364 }
00365
00366 static double local_S41(const JPL_planets planet) {
00367 if (planet == MOON) {
00368 return jpl_file->GetTag("S41M");
00369 } else {
00370 return 0.0;
00371 }
00372 }
00373
00374 static double local_S42(const JPL_planets planet) {
00375 if (planet == MOON) {
00376 return jpl_file->GetTag("S42M");
00377 } else {
00378 return 0.0;
00379 }
00380 }
00381
00382 static double local_S43(const JPL_planets planet) {
00383 if (planet == MOON) {
00384 return jpl_file->GetTag("S43M");
00385 } else {
00386 return 0.0;
00387 }
00388 }
00389
00390 static double local_S44(const JPL_planets planet) {
00391 if (planet == MOON) {
00392 return jpl_file->GetTag("S44M");
00393 } else {
00394 return 0.0;
00395 }
00396 }
00397
00398 JPLBody::JPLBody(const JPL_planets p, const Date & d) : BodyWithEpoch(JPL_planet_name(p),
00399 local_mass(p),
00400 orsa::radius(p),
00401 p,
00402 d,
00403 local_J2(p),
00404 local_J3(p),
00405 local_J4(p),
00406 local_C22(p),
00407 local_C31(p),
00408 local_C32(p),
00409 local_C33(p),
00410 local_C41(p),
00411 local_C42(p),
00412 local_C43(p),
00413 local_C44(p),
00414 local_S31(p),
00415 local_S32(p),
00416 local_S33(p),
00417 local_S41(p),
00418 local_S42(p),
00419 local_S43(p),
00420 local_S44(p)) {
00421 jpl_file->GetEph(epoch,JPLPlanet(),_position,_velocity);
00422 }
00423
00424 JPLBody::JPLBody() : BodyWithEpoch(JPL_planet_name(NONE),local_mass(NONE),orsa::radius(NONE),NONE) {
00425
00426 }
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437 void JPLBody::SetEpoch(const UniverseTypeAwareTime & t) {
00438 epoch = t;
00439 jpl_file->GetEph(epoch,JPLPlanet(),_position,_velocity);
00440 }
00441
00442 void Interpolate(const vector<BodyWithParameter> &b_in, const double x, Body &b_out, Body &err_b_out) {
00443
00444 unsigned int n_points = b_in.size();
00445
00446 Vector p_interpolated, err_p_interpolated;
00447 Vector v_interpolated, err_v_interpolated;
00448
00449 unsigned int j;
00450
00451 vector<VectorWithParameter> pp(n_points), vv(n_points);
00452
00453 for(j=0;j<n_points;j++) {
00454 pp[j].Set(b_in[j].position());
00455 vv[j].Set(b_in[j].velocity());
00456 pp[j].par = vv[j].par = b_in[j].par;
00457 }
00458
00459 Interpolate(pp,x,p_interpolated,err_p_interpolated);
00460 Interpolate(vv,x,v_interpolated,err_v_interpolated);
00461
00462 b_out = b_in[0];
00463 b_out.SetPosition(p_interpolated);
00464 b_out.SetVelocity(v_interpolated);
00465
00466 err_b_out = b_in[0];
00467 err_b_out.SetPosition(err_p_interpolated);
00468 err_b_out.SetVelocity(err_v_interpolated);
00469
00470 }
00471
00472 void print(const Body &b) {
00473 cout << "body name: " << b.name() << " mass: " << b.mass() << endl;
00474
00475
00476
00477
00478
00479 cout << "position: " << b.position().x << " " << b.position().y << " " << b.position().z << endl;
00480 cout << "velocity: " << b.velocity().x << " " << b.velocity().y << " " << b.velocity().z << endl;
00481 }
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598 }