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_file.h"
00026 #include "orsa_error.h"
00027 #include "orsa_secure_math.h"
00028
00029 #include <cstdio>
00030
00031 #include "sdncal.h"
00032 #include "jpleph.h"
00033 #include "jpl_int.h"
00034
00035 #include <iostream>
00036
00037 using namespace std;
00038
00039 namespace orsa {
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 const JPL_planets JPLFile::default_ephem_center = SOLAR_SYSTEM_BARYCENTER;
00050
00051
00052 JPLFile::JPLFile(string name) : calc_velocity(true) {
00053
00054 int N=0;
00055 const int max_N = 256;
00056 char nam[max_N][6];
00057 double val[max_N];
00058
00059
00060 jpl_database = (void *) jpl_init_ephemeris(name.c_str(),NULL,NULL);
00061
00062 if (jpl_database != 0) {
00063 N = ((jpl_eph_data*)jpl_database)->ncon;
00064 if (N > max_N) {
00065 ORSA_ERROR("assumed max_N=%i is smaller than N=%i. Please recompile with a bigger max_N.",max_N,N);
00066 exit(0);
00067 }
00068 jpl_close_ephemeris((jpl_eph_data *)jpl_database);
00069
00070
00071 jpl_database = (void *) jpl_init_ephemeris(name.c_str(),nam,val);
00072 }
00073
00074 if (jpl_database==0) {
00075 ORSA_ERROR("Can't open JPL ephemeris file [%s]",name.c_str());
00076
00077
00078 return;
00079 }
00080
00081
00082
00083
00084
00085 bool_ephem_start_computed = bool_ephem_end_computed = false;
00086
00087 map_tag = new map<std::string,double>;
00088
00089
00090 {
00091 string tag;
00092 char ctag[7];
00093 ctag[6] = 0;
00094
00095 for (int l = 0; l < N; l ++) {
00096 memcpy(ctag, nam[l], 6);
00097 tag = ctag;
00098 remove_leading_trailing_spaces(tag);
00099 (*map_tag)[tag] = val[l];
00100 #if 0
00101 printf(" [l=%03i][%s] = %20.12e\n",l,tag.c_str(),val[l]);
00102 printf(" map_tag[%s] = %20.12e\n",tag.c_str(),(*map_tag)[tag]);
00103 printf(" map_tag[%s] = %20.12e\n",tag.c_str(),GetTag(tag));
00104 #endif
00105 }
00106 }
00107 }
00108
00109 JPLFile::~JPLFile() {
00110 if (jpl_database != 0) jpl_close_ephemeris((jpl_eph_data *)jpl_database);
00111 if (map_tag) delete map_tag;
00112 }
00113
00114 const UniverseTypeAwareTime & JPLFile::EphemStart() {
00115 if (!bool_ephem_start_computed) ComputeEphemStart();
00116 return ephem_start;
00117 }
00118
00119 const UniverseTypeAwareTime & JPLFile::EphemEnd() {
00120 if (!bool_ephem_end_computed) ComputeEphemEnd();
00121 return ephem_end;
00122 }
00123
00124 void JPLFile::ComputeEphemStart() {
00125 jpl_eph_data *jpldb = (jpl_eph_data *) jpl_database;
00126 ephem_start.SetTime(FromUnits(jpldb->ephem_start,DAY));
00127 bool_ephem_start_computed = true;
00128 }
00129
00130 void JPLFile::ComputeEphemEnd() {
00131 jpl_eph_data *jpldb = (jpl_eph_data *) jpl_database;
00132 ephem_end.SetTime(FromUnits(jpldb->ephem_end,DAY));
00133 bool_ephem_end_computed = true;
00134 }
00135
00136 double JPLFile::GetAU_MKS() {
00137 return (GetTag("AU")*1.0e3);
00138 }
00139
00140 double JPLFile::GetMSun_MKS() {
00141
00142 const double au_mks = GetAU_MKS();
00143 const double day_to_second = 24*3600.0;
00144 return (GetTag("GMS")*(au_mks*au_mks*au_mks)/(day_to_second*day_to_second)/GetG_MKS());
00145 }
00146
00147 double JPLFile::GetMJupiter_MKS() {
00148
00149 const double au_mks = GetAU_MKS();
00150 const double day_to_second = 24*3600.0;
00151 return (GetTag("GM5")*(au_mks*au_mks*au_mks)/(day_to_second*day_to_second)/GetG_MKS());
00152 }
00153
00154 double JPLFile::GetMEarth_MKS() {
00155
00156 const double EMRAT = GetTag("EMRAT");
00157
00158
00159
00160 const double au_mks = GetAU_MKS();
00161 const double day_to_second = 24*3600.0;
00162 return ((GetTag("GMB")*EMRAT/(1+EMRAT))*(au_mks*au_mks*au_mks)/(day_to_second*day_to_second)/GetG_MKS());
00163 }
00164
00165 double JPLFile::GetMMoon_MKS() {
00166
00167 const double EMRAT = GetTag("EMRAT");
00168
00169
00170
00171 const double au_mks = GetAU_MKS();
00172 const double day_to_second = 24*3600.0;
00173 return ((GetTag("GMB")/(1+EMRAT))*(au_mks*au_mks*au_mks)/(day_to_second*day_to_second)/GetG_MKS());
00174 }
00175
00176 double JPLFile::GetC_MKS() {
00177 return (GetTag("CLIGHT")*1.0e3);
00178 }
00179
00180 double JPLFile::GetREarth_MKS() {
00181 return (GetTag("RE")*1.0e3);
00182 }
00183
00184 double JPLFile::GetRMoon_MKS() {
00185 return (GetTag("AM")*1.0e3);
00186 }
00187
00188 double JPLFile::GetMass(const JPL_planets planet) {
00189
00190
00191
00192
00193
00194 const double EMRAT = GetTag("EMRAT");
00195
00196 double GM=0;
00197
00198 switch(planet) {
00199 case MERCURY:
00200 GM = GetTag("GM1");
00201 break;
00202 case VENUS:
00203 GM = GetTag("GM2");
00204 break;
00205 case MARS:
00206 GM = GetTag("GM4");
00207 break;
00208 case JUPITER:
00209 GM = GetTag("GM5");
00210 break;
00211 case SATURN:
00212 GM = GetTag("GM6");
00213 break;
00214 case URANUS:
00215 GM = GetTag("GM7");
00216 break;
00217 case NEPTUNE:
00218 GM = GetTag("GM8");
00219 break;
00220 case PLUTO:
00221 GM = GetTag("GM9");
00222 break;
00223 case EARTH:
00224 GM = GetTag("GMB")*EMRAT/(1+EMRAT);
00225 break;
00226 case MOON:
00227 GM = GetTag("GMB")/(1+EMRAT);
00228 break;
00229 case EARTH_MOON_BARYCENTER:
00230 GM = GetTag("GMB");
00231 break;
00232 case SUN:
00233 GM = GetTag("GMS");
00234 break;
00235 default:
00236 GM = 0;
00237 break;
00238 }
00239
00240
00241 GM = FromUnits(FromUnits(GM,AU,3),DAY,-2);
00242
00243
00244 return (GM/GetG());
00245 }
00246
00247 double JPLFile::GetTag(string tag) {
00248 remove_leading_trailing_spaces(tag);
00249 return (*map_tag)[tag];
00250 }
00251
00252
00253
00254 void JPLFile::GetEph(const UniverseTypeAwareTime & date, JPL_planets target, JPL_planets center, Vector & position, Vector & velocity) {
00255 double xv[6];
00256
00257 jpl_eph_data * jpldb = (jpl_eph_data *) jpl_database;
00258
00259 if ( (date < EphemStart()) ||
00260 (date > EphemEnd()) ) {
00261 ORSA_WARNING("requested time out of the jpl database range");
00262 return;
00263 }
00264
00265 jpl_pleph(jpldb,date.GetDate().GetJulian(ET),target,center,xv,calc_velocity ? 1 : 0);
00266
00267 if ((target==NUTATIONS) ||
00268 (target==LIBRATIONS)) {
00269
00270 position.Set(xv[0],xv[1],xv[2]);
00271 velocity.Set(xv[3],xv[4],xv[5]);
00272 return;
00273 }
00274
00275
00276 xv[0] = FromUnits(xv[0],AU);
00277 xv[1] = FromUnits(xv[1],AU);
00278 xv[2] = FromUnits(xv[2],AU);
00279
00280 position.Set(xv[0],xv[1],xv[2]);
00281
00282 if (calc_velocity) {
00283
00284 xv[3] = FromUnits(xv[3],AU);
00285 xv[4] = FromUnits(xv[4],AU);
00286 xv[5] = FromUnits(xv[5],AU);
00287
00288 xv[3] = FromUnits(xv[3],DAY,-1);
00289 xv[4] = FromUnits(xv[4],DAY,-1);
00290 xv[5] = FromUnits(xv[5],DAY,-1);
00291
00292 velocity.Set(xv[3],xv[4],xv[5]);
00293 }
00294
00295
00296
00297 if (universe->GetReferenceSystem() == ECLIPTIC) {
00298 Angle obl = obleq_J2000();
00299 position.rotate(0.0,-obl.GetRad(),0.0);
00300 velocity.rotate(0.0,-obl.GetRad(),0.0);
00301 }
00302 }
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413 void SetupSolarSystem(Frame & frame, const list<JPL_planets> & l, const UniverseTypeAwareTime & t) {
00414
00415 frame.clear();
00416 frame.SetTime(t);
00417
00418
00419
00420
00421
00422
00423
00424
00425 if (t < jpl_file->EphemStart()) {
00426 ORSA_WARNING("epoch requested is before ephem file start time! (%g < %g)",t.GetTime(),jpl_file->EphemStart().GetTime());
00427 return;
00428 }
00429
00430 if (t > jpl_file->EphemEnd()) {
00431 ORSA_WARNING("epoch requested is after ephem file end time! (%g > %g)",t.GetTime(),jpl_file->EphemStart().GetTime());
00432 return;
00433 }
00434
00435
00436
00437 frame.push_back(jpl_cache->GetJPLBody(SUN,t));
00438
00439 JPL_planets pl;
00440 list<JPL_planets>::const_iterator it = l.begin();
00441 while (it != l.end()) {
00442 pl = *it;
00443 if (pl==SUN) {
00444 ++it;
00445 continue;
00446 } else if (pl==EARTH_AND_MOON) {
00447
00448 frame.push_back(jpl_cache->GetJPLBody(EARTH,t));
00449
00450 frame.push_back(jpl_cache->GetJPLBody(MOON,t));
00451 } else {
00452
00453 frame.push_back(jpl_cache->GetJPLBody(pl,t));
00454 }
00455 ++it;
00456 }
00457 }
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480 string JPL_planet_name(const JPL_planets p) {
00481
00482 string name;
00483
00484 switch(p) {
00485 case SUN: name = "Sun"; break;
00486 case MERCURY: name = "Mercury"; break;
00487 case VENUS: name = "Venus"; break;
00488 case EARTH: name = "Earth"; break;
00489 case MARS: name = "Mars"; break;
00490 case JUPITER: name = "Jupiter"; break;
00491 case SATURN: name = "Saturn"; break;
00492 case URANUS: name = "Uranus"; break;
00493 case NEPTUNE: name = "Neptune"; break;
00494 case PLUTO: name = "Pluto"; break;
00495
00496 case MOON: name = "Moon"; break;
00497 case EARTH_MOON_BARYCENTER: name = "Earth-Moon barycenter"; break;
00498
00499 default: break;
00500 }
00501
00502 return (name);
00503 }
00504
00505
00506
00507 double radius(const JPL_planets p) {
00508 double r=0;
00509 switch(p) {
00510 case SUN: r = FromUnits(695000. ,KM); break;
00511 case MERCURY: r = FromUnits( 2440. ,KM); break;
00512 case VENUS: r = FromUnits( 6051.84,KM); break;
00513 case EARTH: r = FromUnits( 6371.01,KM); break;
00514 case MARS: r = FromUnits( 3389.92,KM); break;
00515 case JUPITER: r = FromUnits( 69911. ,KM); break;
00516 case SATURN: r = FromUnits( 58232. ,KM); break;
00517 case URANUS: r = FromUnits( 25362. ,KM); break;
00518 case NEPTUNE: r = FromUnits( 24624. ,KM); break;
00519 case PLUTO: r = FromUnits( 1151. ,KM); break;
00520
00521 case MOON: r = FromUnits( 1738. ,KM); break;
00522
00523 default: break;
00524 }
00525 return (r);
00526 }
00527
00528
00529
00530 JPLCache::JPLCache() {
00531
00532 }
00533
00534 JPLCache::~JPLCache() {
00535
00536 }
00537
00538 const JPLBody & JPLCache::GetJPLBody(const JPL_planets p, const UniverseTypeAwareTime & t) {
00539 data_map_type & data = big_map[p];
00540 data_map_type::const_iterator it = data.find(t);
00541 if (it != data.end()) {
00542
00543 return ((*it).second);
00544 } else {
00545
00546
00547
00548 return (data[t] = JPLBody(p,t));
00549 }
00550 }
00551
00552 void JPLCache::Clear() {
00553 big_map.clear();
00554 }
00555
00556 }