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