orsa_file_jpl.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_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   // UPDATE (Apr 2004): the results are the same, no matter what default object is selected
00043   // const JPL_planets JPLFile::default_ephem_center = SOLAR_SYSTEM_BARYCENTER;
00044   // const JPL_planets JPLFile::default_ephem_center = SUN;
00045   //
00046   // CORRECT: keep the EARTH as default center
00047   // const JPL_planets JPLFile::default_ephem_center = EARTH;
00048   //
00049   const JPL_planets JPLFile::default_ephem_center = SOLAR_SYSTEM_BARYCENTER;
00050   
00051   // JPL
00052   JPLFile::JPLFile(string name) : calc_velocity(true) {
00053     
00054     int N=0; // number of elements in nam[][6] and val[]
00055     const int max_N = 256;
00056     char nam[max_N][6];
00057     double val[max_N];
00058     
00059     // jpl_database = (jpl_eph_data *) jpl_init_ephemeris(name.c_str(),nam,val);
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       // jpl_database = (void *) jpl_init_ephemeris(name.c_str(),nam,val);
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       // exit is too much, should handle better this problem
00077       // exit(0);
00078       return;
00079     }
00080     
00081     // jpl_eph_data *jpldb = (jpl_eph_data *) jpl_database;
00082     // ephem_start.SetTime(FromUnits(jpldb->ephem_start,DAY));
00083     // ephem_end.SetTime(  FromUnits(jpldb->ephem_end,  DAY));
00084     
00085     bool_ephem_start_computed = bool_ephem_end_computed = false;  
00086     
00087     map_tag = new map<std::string,double>;
00088     
00089     // fill map_tag, DON'T REMOVE
00090     {
00091       string tag;
00092       char ctag[7];
00093       ctag[6] = 0;
00094       // cerr << "N: " << N << endl;
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     //return (GetTag("GMS")*secure_pow(GetAU_MKS(),3.0)*secure_pow(24*3600.0,-2.0)/GetG_MKS());
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     // return (GetTag("GM5")*secure_pow(GetAU_MKS(),3.0)*secure_pow(24*3600.0,-2.0)/GetG_MKS());
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     // const double EMRAT = val[7]; // Earth/Moon mass ratio
00156     const double EMRAT = GetTag("EMRAT");
00157     // return ((val[10]*EMRAT/(1+EMRAT))*pow(AU_MKS,3)*pow(24*3600.0,-2)/G_MKS);
00158     // return ((GetTag("GMB")*EMRAT/(1+EMRAT))*pow(AU_MKS,3.0)*pow(24*3600.0,-2.0)/G_MKS);
00159     // return ((GetTag("GMB")*EMRAT/(1+EMRAT))*secure_pow(GetAU_MKS(),3.0)*secure_pow(24*3600.0,-2.0)/GetG_MKS());
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     // const double EMRAT = val[7]; // Earth/Moon mass ratio
00167     const double EMRAT = GetTag("EMRAT");
00168     // return ((val[10]/(1+EMRAT))*pow(AU_MKS,3)*pow(24*3600.0,-2)/G_MKS);
00169     // return ((GetTag("GMB")/(1+EMRAT))*pow(AU_MKS,3.0)*pow(24*3600.0,-2.0)/G_MKS);
00170     // return ((GetTag("GMB")/(1+EMRAT))*secure_pow(GetAU_MKS(),3.0)*secure_pow(24*3600.0,-2.0)/GetG_MKS());
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     // const double file_G = FromUnits( FromUnits( FromUnits(val[17],AU,3),MSUN,-1),DAY,-2);
00191     // cerr << "JPLFile::GetMass G=" << file_G << endl;
00192     
00193     // const double EMRAT = val[7]; // Earth/Moon mass ratio
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     // take the right units
00241     GM = FromUnits(FromUnits(GM,AU,3),DAY,-2);
00242     
00243     // return (GM/file_G);
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       // no units correction needed, are radians
00270       position.Set(xv[0],xv[1],xv[2]);
00271       velocity.Set(xv[3],xv[4],xv[5]);
00272       return;
00273     }
00274     
00275     // position (from AU)
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       // velocity (from AU/DAY)
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     // THIS IS THE ONLY CORRECT ROTATION:
00296     // from mean equatorial J2000 (ICRF) to mean ecliptic J2000
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   // small utility
00305   /* 
00306      static void SetMRV(double &mass, Vector &r, Vector &v, JPLFile &jf, Date d, JPL_planets p) {
00307      mass = jf.GetMass(p);  
00308      jf.GetEph(d,p,r,v);
00309      }
00310   */
00311   
00312   /* 
00313      static void SetMRV(double &mass, Vector &r, Vector &v, JPLFile * jf, const UniverseTypeAwareTime & t, JPL_planets p) {
00314      mass = jf->GetMass(p);  
00315      jf->GetEph(t,p,r,v);
00316      }
00317   */
00318   
00319   //! The Sun is automatically included
00320   /* 
00321      void SetupSolarSystem(Frame & frame, const list<JPL_planets> & l, const UniverseTypeAwareTime & t) {
00322      
00323      // fprintf(stderr,"sss: ls=%i t=%f\n",l.size(),t.GetTime());
00324      
00325      const Date epoch(t.GetDate());
00326      
00327      frame.clear();
00328      frame.SetDate(epoch);
00329      
00330      const string jpl_path = config->paths[JPL_EPHEM_FILE]->GetValue().c_str();
00331      JPLFile jf(jpl_path);
00332      
00333      string name;
00334      double mass;
00335      Vector r,v;
00336      
00337      // date checks
00338      
00339      if (t.GetTime() < jf.EphemStart().GetTime()) { 
00340      ORSA_WARNING("epoch requested is before ephem file start time! (%g < %g)",t.GetTime(),jf.EphemStart().GetTime());
00341      return; 
00342      }
00343      //
00344      if (t.GetTime() > jf.EphemEnd().GetTime()) {
00345      ORSA_WARNING("epoch requested is after ephem file end time! (%g > %g)",t.GetTime(),jf.EphemStart().GetTime());
00346      return;
00347      }
00348      
00349      // auto include the Sun
00350      name="Sun"; SetMRV(mass,r,v,jf,epoch,SUN); frame.push_back(Body(name,mass,r,v));
00351      
00352      JPL_planets pl;
00353      list<JPL_planets>::const_iterator it = l.begin();
00354      while (it != l.end()) {
00355      pl = *it;
00356      if (pl==SUN) {
00357      ++it;
00358      continue; // sun already included
00359      } else if (pl==EARTH_AND_MOON) {
00360      pl = EARTH; name = JPL_planet_name(pl); SetMRV(mass,r,v,jf,epoch,pl);  frame.push_back(Body(name,mass,r,v));
00361      pl = MOON;  name = JPL_planet_name(pl); SetMRV(mass,r,v,jf,epoch,pl);  frame.push_back(Body(name,mass,r,v));
00362      } else {
00363      name = JPL_planet_name(pl); SetMRV(mass,r,v,jf,epoch,pl); frame.push_back(Body(name,mass,r,v));
00364      }
00365      ++it;
00366      }    
00367      }
00368   */
00369   
00370   //! The Sun is automatically included
00371   /* 
00372      void SetupSolarSystem(Frame & frame, const list<JPL_planets> & l, const UniverseTypeAwareTime & t) {
00373      
00374      frame.clear();
00375      frame.SetTime(t);
00376      
00377      string name;
00378      double mass;
00379      Vector r,v;
00380      
00381      // date checks
00382      if (t < jpl_file->EphemStart()) { 
00383      ORSA_WARNING("epoch requested is before ephem file start time! (%g < %g)",t.GetTime(),jpl_file->EphemStart().GetTime());
00384      return; 
00385      }
00386      //
00387      if (t > jpl_file->EphemEnd()) {
00388      ORSA_WARNING("epoch requested is after ephem file end time! (%g > %g)",t.GetTime(),jpl_file->EphemStart().GetTime());
00389      return;
00390      }
00391      
00392      // auto include the Sun
00393      name="Sun"; SetMRV(mass,r,v,jpl_file,t,SUN); frame.push_back(Body(name,mass,r,v));
00394      
00395      JPL_planets pl;
00396      list<JPL_planets>::const_iterator it = l.begin();
00397      while (it != l.end()) {
00398      pl = *it;
00399      if (pl==SUN) {
00400      ++it;
00401      continue; // sun already included
00402      } else if (pl==EARTH_AND_MOON) {
00403      pl = EARTH; name = JPL_planet_name(pl); SetMRV(mass,r,v,jpl_file,t,pl);  frame.push_back(Body(name,mass,r,v));
00404      pl = MOON;  name = JPL_planet_name(pl); SetMRV(mass,r,v,jpl_file,t,pl);  frame.push_back(Body(name,mass,r,v));
00405      } else {
00406      name = JPL_planet_name(pl); SetMRV(mass,r,v,jpl_file,t,pl); frame.push_back(Body(name,mass,r,v));
00407      }
00408      ++it;
00409      }    
00410      }
00411   */
00412   //! The Sun is automatically included
00413   void SetupSolarSystem(Frame & frame, const list<JPL_planets> & l, const UniverseTypeAwareTime & t) {
00414     
00415     frame.clear();
00416     frame.SetTime(t);
00417     
00418     /* 
00419        string name;
00420        double mass;
00421        Vector r,v;
00422     */
00423     
00424     // date checks
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     // auto include the Sun
00436     // name="Sun"; SetMRV(mass,r,v,jpl_file,t,SUN); frame.push_back(Body(name,mass,r,v));
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; // sun already included
00446       } else if (pl==EARTH_AND_MOON) {
00447         // pl = EARTH; name = JPL_planet_name(pl); SetMRV(mass,r,v,jpl_file,t,pl);  frame.push_back(Body(name,mass,r,v));
00448         frame.push_back(jpl_cache->GetJPLBody(EARTH,t));
00449         // pl = MOON;  name = JPL_planet_name(pl); SetMRV(mass,r,v,jpl_file,t,pl);  frame.push_back(Body(name,mass,r,v));
00450         frame.push_back(jpl_cache->GetJPLBody(MOON,t));
00451       } else {
00452         // name = JPL_planet_name(pl); SetMRV(mass,r,v,jpl_file,t,pl); frame.push_back(Body(name,mass,r,v));
00453         frame.push_back(jpl_cache->GetJPLBody(pl,t));
00454       }
00455       ++it;
00456     }    
00457   }
00458   
00459   ////
00460   
00461   /* 
00462      bool is_genuine_planet(const JPL_planets p) {
00463      bool is_genuine;
00464      switch(p) {
00465      case SOLAR_SYSTEM_BARYCENTER:
00466      case EARTH_MOON_BARYCENTER:
00467      case NUTATIONS:
00468      case LIBRATIONS:
00469      case EARTH_AND_MOON:
00470      is_genuine = false;
00471      break;
00472      default:
00473      is_genuine = true;
00474      break;
00475      }
00476      return (is_genuine);
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   // data from http://ssd.jpl.nasa.gov/phys_props_planets.html (Dec 6 2003)
00506   // Sun data: http://www.solarviews.com/eng/sun.htm
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   // JPLCache 
00529   
00530   JPLCache::JPLCache() {
00531     // jf = new JPLFile(config->paths[JPL_EPHEM_FILE]->GetValue().c_str());
00532   }
00533   
00534   JPLCache::~JPLCache() {
00535     // if (jf) delete jf;
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       // ORSA_ERROR("JPLCache::GetJPLBody(...) ==> Found something in cache...");
00543       return ((*it).second);
00544     } else {
00545       // ORSA_ERROR("JPLCache::GetJPLBody(...) ==> Adding object to cache...");
00546       // data[t] = JPLBody(p,t);
00547       // return data[t];
00548       return (data[t] = JPLBody(p,t));
00549     }
00550   }
00551   
00552   void JPLCache::Clear() {
00553     big_map.clear();
00554   }
00555   
00556 } // namespace orsa

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