orsa_units.h

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 #ifndef _ORSA_UNITS_H
00026 #define _ORSA_UNITS_H
00027 
00028 #include <cmath>
00029 #include <string>
00030 #include <cstdio>
00031 #include <cstdlib>
00032 
00033 #include "orsa_secure_math.h"
00034 #include "orsa_coord.h"
00035 #include "orsa_error.h"
00036 
00037 namespace orsa {
00038   
00039   enum time_unit {
00040     YEAR=1,
00041     DAY=2,
00042     HOUR=3,
00043     MINUTE=4,
00044     SECOND=5
00045   };
00046   
00047   inline void convert(time_unit &tu, const unsigned int i)  {
00048     switch(i) {
00049     case 1: tu = YEAR;   break;
00050     case 2: tu = DAY;    break;
00051     case 3: tu = HOUR;   break;
00052     case 4: tu = MINUTE; break;
00053     case 5: tu = SECOND; break;
00054       //
00055     default:
00056       ORSA_ERROR("conversion problem: i = %i",i);
00057       break;       
00058     }
00059   }
00060   
00061   enum length_unit {
00062     MPARSEC=1,
00063     KPARSEC=2,
00064     PARSEC=3,
00065     LY=4,
00066     AU=5,
00067     EARTHMOON=6,
00068     REARTH=7,
00069     RMOON=8,
00070     KM=9,
00071     M=10,
00072     CM=11,
00073     // aliases
00074     LD=EARTHMOON,
00075     ER=REARTH,
00076     MR=RMOON
00077   };
00078   
00079   inline void convert(length_unit &lu, const unsigned int i)  {
00080     switch(i) {
00081     case 1:  lu = MPARSEC;   break;
00082     case 2:  lu = KPARSEC;   break;
00083     case 3:  lu = PARSEC;    break;
00084     case 4:  lu = LY;        break;
00085     case 5:  lu = AU;        break;
00086     case 6:  lu = EARTHMOON; break;
00087     case 7:  lu = REARTH;    break;
00088     case 8:  lu = RMOON;     break;
00089     case 9:  lu = KM;        break;
00090     case 10: lu = M;         break;
00091     case 11: lu = CM;        break;
00092       //
00093     default:
00094       ORSA_ERROR("conversion problem: i = %i",i);
00095       break;       
00096     }
00097   }
00098   
00099   enum mass_unit {
00100     MSUN=1,
00101     MJUPITER=2,
00102     MEARTH=3,
00103     MMOON=4,
00104     KG=5,
00105     GRAM=6
00106   };
00107   
00108   inline void convert(mass_unit &mu, const unsigned int i)  {
00109     switch(i) {
00110     case 1: mu = MSUN;     break;
00111     case 2: mu = MJUPITER; break;
00112     case 3: mu = MEARTH;   break;
00113     case 4: mu = MMOON;    break;
00114     case 5: mu = KG;       break;
00115     case 6: mu = GRAM;     break;    
00116       //
00117     default:
00118       ORSA_ERROR("conversion problem: i = %i",i);
00119       break;       
00120     }
00121   }
00122   
00123   template <class UNIT> class UnitBaseScale {
00124   public:
00125     UnitBaseScale() {}
00126     UnitBaseScale(UNIT unit) { base_unit = unit; }
00127  
00128   public:
00129     void Set(UNIT unit) { base_unit = unit; }
00130     const UNIT & GetBaseUnit() const { return base_unit; }
00131     
00132   private:
00133     UNIT base_unit;
00134   };
00135   
00136   class Units {
00137     
00138   private:
00139     UnitBaseScale<time_unit>   Time;
00140     UnitBaseScale<length_unit> Length;
00141     UnitBaseScale<mass_unit>   Mass;
00142     
00143   public:
00144     Units();
00145     Units(time_unit, length_unit, mass_unit);
00146     
00147   private:
00148     void init_base();
00149     
00150   private:
00151     void TryToSetUnitsFromJPLFile();
00152     
00153   public:
00154     void SetSystem(time_unit tu, length_unit lu, mass_unit mu);
00155     
00156     // constants
00157     inline double GetG()    const { return G; };
00158     inline double GetMSun() const { return MSun; };
00159     inline double GetC()    const { return c; }; // c = speed of light
00160     
00161     double GetG_MKS() const;
00162     
00163   public:
00164     std::string label(time_unit)   const;
00165     std::string label(length_unit) const;
00166     std::string label(mass_unit)   const;
00167     
00168   public:
00169     inline std::string TimeLabel()   const { return label(  Time.GetBaseUnit()); };
00170     inline std::string LengthLabel() const { return label(Length.GetBaseUnit()); };
00171     inline std::string MassLabel()   const { return label(  Mass.GetBaseUnit()); };
00172     
00173     // conversions
00174   public:
00175     /* 
00176        inline double FromUnits(double x, time_unit   t_in, double power=1.0) const { return (x*secure_pow(GetTimeScale(t_in)/GetTimeScale(),power));     }
00177        inline double FromUnits(double x, length_unit l_in, double power=1.0) const { return (x*secure_pow(GetLengthScale(l_in)/GetLengthScale(),power)); }
00178        inline double FromUnits(double x, mass_unit   m_in, double power=1.0) const { return (x*secure_pow(GetMassScale(m_in)/GetMassScale(),power));     }
00179     */
00180     
00181   private:
00182     inline static double __int_pow__(const double x, const int p) {
00183       if (p == 0) return 1.0;
00184       double _pow = x;
00185       const unsigned int max_k = static_cast<unsigned int>(std::abs(p));
00186       for (unsigned int k=1; k < max_k; ++k) {
00187         _pow *= x;
00188       }
00189       if (p < 0) _pow = 1.0/_pow;
00190       return _pow;
00191     }
00192     
00193     // conversions
00194   public:
00195     inline double FromUnits(const double x, const time_unit   t_in, const int power = 1) const { return (x*__int_pow__(GetTimeScale(t_in)/GetTimeScale(),power));     }
00196     inline double FromUnits(const double x, const length_unit l_in, const int power = 1) const { return (x*__int_pow__(GetLengthScale(l_in)/GetLengthScale(),power)); }
00197     inline double FromUnits(const double x, const mass_unit   m_in, const int power = 1) const { return (x*__int_pow__(GetMassScale(m_in)/GetMassScale(),power));     }
00198     
00199   public:
00200     inline const time_unit   & GetTimeBaseUnit()   const { return Time.GetBaseUnit(); };
00201     inline const length_unit & GetLengthBaseUnit() const { return Length.GetBaseUnit(); };
00202     inline const mass_unit   & GetMassBaseUnit()   const { return Mass.GetBaseUnit(); };
00203     
00204   protected:
00205     double GetTimeScale(  const   time_unit tu) const; 
00206     double GetLengthScale(const length_unit lu) const; 
00207     double GetMassScale(  const   mass_unit mu) const; 
00208     
00209     double GetTimeScale()   const;
00210     double GetLengthScale() const; 
00211     double GetMassScale()   const; 
00212     
00213     void Recompute();
00214     
00215   private:
00216     double G,G_base;
00217     double MSun,MSun_base;
00218     double MJupiter_base, MEarth_base, MMoon_base;
00219     double AU_base;
00220     double c,c_base;
00221     double r_earth_base;
00222     double r_moon_base;
00223     double parsec_base;
00224   };
00225   
00226   // world visible units // defined in orsa_universe.cc
00227   extern Units * units;
00228   
00229   // defined in orsa_universe.cc
00230   double GetG();
00231   double GetG_MKS();
00232   double GetMSun();
00233   double GetC();
00234   
00235   // everything in orsa_universe.cc
00236   /* 
00237      double FromUnits(double, time_unit,   double = 1.0);
00238      double FromUnits(double, length_unit, double = 1.0);
00239      double FromUnits(double, mass_unit,   double = 1.0); 
00240   */
00241   //
00242   double FromUnits(const double, const time_unit,   const int = 1);
00243   double FromUnits(const double, const length_unit, const int = 1);
00244   double FromUnits(const double, const mass_unit,   const int = 1); 
00245   //
00246   std::string TimeLabel();
00247   std::string LengthLabel();
00248   std::string MassLabel();
00249   
00250   //! TimeScale enum, useful only when using a Real Universe. More information can be obtained here: http://www.hartrao.ac.za/nccsdoc/slalib/sun67.htx/node217.html
00251   enum TimeScale {
00252     UTC=1,
00253     UT=2,
00254     TAI=3,
00255     TDT=4,
00256     GPS=5,
00257     // aliases
00258     UT1=UT,
00259     ET=TDT,
00260     TT=TDT
00261   };
00262   
00263   inline void convert(TimeScale &ts, const unsigned int i)  {
00264     switch(i) {
00265     case 1: ts = UTC; break;
00266     case 2: ts = UT;  break;
00267     case 3: ts = TAI; break;
00268     case 4: ts = TDT; break;
00269     case 5: ts = GPS; break;
00270       //
00271     default:
00272       ORSA_ERROR("conversion problem: i = %i",i);
00273       break;       
00274     }
00275   }
00276   
00277   extern TimeScale default_Date_timescale;
00278   
00279   std::string TimeScaleLabel(TimeScale);
00280   
00281   class UniverseTypeAwareTime;
00282   
00283   class UniverseTypeAwareTimeStep;
00284   
00285   class TimeStep {
00286   public:
00287     TimeStep();
00288     TimeStep(const unsigned int days, const unsigned int day_fraction, const int sign);
00289     TimeStep(const double t);
00290     TimeStep(const TimeStep &);
00291     
00292   public:
00293     double GetDouble() const;
00294     
00295   public:
00296     TimeStep & operator += (const TimeStep &);
00297     TimeStep & operator -= (const TimeStep &);
00298     
00299     TimeStep operator + (const TimeStep &) const;
00300     TimeStep operator - (const TimeStep &) const;
00301     
00302     TimeStep & operator *= (const int);
00303     TimeStep & operator *= (const double);
00304     
00305   public:
00306     void AddDays(const unsigned int, const int);
00307     void AddDayFractions(const unsigned int, const int);
00308     
00309     // sign operators 
00310     TimeStep operator + () const;
00311     TimeStep operator - () const;
00312     
00313     bool operator == (const TimeStep &) const;
00314     bool operator != (const TimeStep &) const;
00315     
00316     bool operator > (const TimeStep &) const;
00317     bool operator < (const TimeStep &) const;
00318     
00319     inline bool operator >= (const TimeStep & ts) const {
00320       if ((*this) == ts) return true;
00321       if ((*this) >  ts) return true;
00322       return false;
00323     }
00324     
00325     inline bool operator <= (const TimeStep & ts) const {
00326       if ((*this) == ts) return true;
00327       if ((*this) <  ts) return true;
00328       return false;      
00329     }
00330     
00331     TimeStep absolute() const;
00332     
00333     bool IsZero() const;
00334     
00335     static unsigned int max_day_fraction() {
00336       return 864000000; // 86400 ( = second per day) * 10000 ( = 1.0e-4 seconds resolution)
00337     }
00338     
00339   private:
00340     void internal_check();
00341     
00342   public:
00343     unsigned int days()         const { return _days;         } 
00344     unsigned int day_fraction() const { return _day_fraction; }
00345     int          sign()         const { return _sign;         }
00346     
00347   private:
00348     unsigned int _days;         // always positive
00349     unsigned int _day_fraction; // always positive
00350     int          _sign;         // +1 or -1
00351   };
00352   
00353   /*! A note on Gregorian and Julian calendars:
00354    * since its introduction in 1582, the Gregorian calendar
00355    * has been adopted in different years from different countries,
00356    * so the date obtained from the Date class can be different from 
00357    * the real 'old' one which used the Julian date, and usually
00358    * the difference is of a few days.
00359    *
00360    * In particular, ORSA applies the Gregorian calendar in any epoch,
00361    * i.e. even before 1582. This appears to be the simplest solution,
00362    * at least for the moment.
00363    * 
00364    * For more info, check out i.e. http://www.dome-igm.com/convers.htm
00365    */
00366   
00367   class Date {
00368   public:
00369     Date();
00370     Date(const Date &);
00371     Date(const UniverseTypeAwareTime &);
00372     
00373   public:
00374     void SetGregor(int   y, int   m, double  d,                              TimeScale ts = default_Date_timescale);
00375     void SetGregor(int   y, int   m, int     d, int H, int M, int S, int ms, TimeScale ts = default_Date_timescale);
00376     //
00377     void GetGregor(int & y, int & m, int &   d, TimeScale ts = default_Date_timescale) const;
00378     double GetDayFraction(                      TimeScale ts = default_Date_timescale) const;
00379     unsigned int GetDayFraction_unsigned_int(   TimeScale ts = default_Date_timescale) const;
00380     //
00381     void   SetJulian(double,               TimeScale ts = default_Date_timescale);
00382     void   GetJulian(double &,             TimeScale ts = default_Date_timescale) const;
00383     double GetJulian(                      TimeScale ts = default_Date_timescale) const;
00384     // 
00385     double Time()    const;
00386     double GetTime() const;
00387     void   SetTime(const double);
00388     
00389     Date & operator += (const UniverseTypeAwareTimeStep &);
00390     Date & operator -= (const UniverseTypeAwareTimeStep &);
00391     
00392     bool operator == (const Date &) const;
00393     bool operator != (const Date &) const;
00394     
00395     bool operator < (const Date &) const;
00396     bool operator > (const Date &) const;
00397     
00398   public:
00399     void SetNow();
00400     void SetToday();
00401     void SetJ2000();
00402     
00403   public:
00404     //! delta_seconds = to - from   ('to' minus 'from')
00405     static double delta_seconds(int y, int m, int d, const TimeScale from, const TimeScale to=default_Date_timescale);
00406     
00407     //! returns a TimeStep as big as the Date
00408     TimeStep GetTimeStep() const {
00409       return TimeStep(sdn,df,+1);
00410     }
00411     
00412   private:
00413     unsigned int sdn;  // internal day format
00414     unsigned int  df;  // day fraction; compare it with TimeStep::max_day_fraction()
00415   };
00416   
00417   /* 
00418      inline bool operator == (const Date &x, const Date &y) {
00419      return (x.GetTimeStep() == y.GetTimeStep());
00420      }
00421      
00422      inline bool operator!= (const Date &x, const Date &y) {
00423      return !(x==y);
00424      }
00425   */
00426   
00427   // UniverseTypeAwareTime
00428   
00429   //! A Date like class which represents time
00430   //! as Date in Real Universes and as a double
00431   //! in Simulated Universes  
00432   class UniverseTypeAwareTime {
00433   public:
00434     UniverseTypeAwareTime();
00435     UniverseTypeAwareTime(const double);
00436     UniverseTypeAwareTime(const Date&);
00437     UniverseTypeAwareTime(const UniverseTypeAwareTime&);
00438     
00439   public:
00440     double Time()    const;
00441     double GetTime() const;
00442     Date   GetDate() const;
00443     virtual void SetTime(const double);
00444     virtual void SetTime(const Date &);
00445     virtual void SetDate(const Date &);
00446     virtual void SetTime(const UniverseTypeAwareTime &);
00447     
00448     UniverseTypeAwareTime & operator += (const UniverseTypeAwareTimeStep &);
00449     UniverseTypeAwareTime & operator -= (const UniverseTypeAwareTimeStep &);
00450     
00451     // ***
00452     // UniverseTypeAwareTimeStep operator + (const UniverseTypeAwareTime &) const;
00453     UniverseTypeAwareTimeStep operator - (const UniverseTypeAwareTime &) const;
00454     
00455     // UniverseTypeAwareTimeStep operator + (const UniverseTypeAwareTimeStep &) const;
00456     // UniverseTypeAwareTimeStep operator - (const UniverseTypeAwareTimeStep &) const;
00457     
00458     UniverseTypeAwareTime operator + (const UniverseTypeAwareTimeStep &) const;
00459     UniverseTypeAwareTime operator - (const UniverseTypeAwareTimeStep &) const;
00460     
00461     bool operator == (const UniverseTypeAwareTime &) const;
00462     bool operator != (const UniverseTypeAwareTime &) const;
00463     
00464     bool operator < (const UniverseTypeAwareTime &) const;
00465     bool operator > (const UniverseTypeAwareTime &) const;
00466     
00467     bool operator <= (const UniverseTypeAwareTime &) const;
00468     bool operator >= (const UniverseTypeAwareTime &) const;
00469     
00470   protected:
00471     double time;
00472     Date   date;
00473   };
00474   
00475   class UniverseTypeAwareTimeStep {
00476   public:
00477     UniverseTypeAwareTimeStep();
00478     UniverseTypeAwareTimeStep(const double);
00479     UniverseTypeAwareTimeStep(const TimeStep &);
00480     UniverseTypeAwareTimeStep(const int days, const unsigned int day_fraction, const int sign = +1);
00481     UniverseTypeAwareTimeStep(const UniverseTypeAwareTimeStep &);
00482     
00483   public:
00484     UniverseTypeAwareTimeStep & operator += (const UniverseTypeAwareTimeStep &);
00485     UniverseTypeAwareTimeStep & operator -= (const UniverseTypeAwareTimeStep &);
00486     UniverseTypeAwareTimeStep & operator *= (const int);
00487     UniverseTypeAwareTimeStep & operator *= (const double);
00488     
00489   public:
00490     double GetDouble() const;
00491     
00492   public:
00493     UniverseTypeAwareTimeStep operator + () const;
00494     UniverseTypeAwareTimeStep operator - () const;
00495   
00496   public:
00497     UniverseTypeAwareTimeStep operator + (const UniverseTypeAwareTimeStep &) const;
00498     UniverseTypeAwareTimeStep operator - (const UniverseTypeAwareTimeStep &) const;
00499     
00500     UniverseTypeAwareTimeStep operator + (const UniverseTypeAwareTime &) const;
00501     UniverseTypeAwareTimeStep operator - (const UniverseTypeAwareTime &) const;
00502     
00503     UniverseTypeAwareTimeStep absolute() const;
00504     
00505     bool IsZero() const;
00506     
00507     bool operator < (const UniverseTypeAwareTimeStep &) const;
00508     bool operator > (const UniverseTypeAwareTimeStep &) const;
00509     
00510     bool operator < (const UniverseTypeAwareTime &) const;
00511     bool operator > (const UniverseTypeAwareTime &) const;
00512     
00513     bool operator < (const double) const;
00514     bool operator > (const double) const;
00515     
00516     bool operator == (const UniverseTypeAwareTimeStep &) const;
00517     bool operator != (const UniverseTypeAwareTimeStep &) const;
00518     
00519   public:
00520     unsigned int days()         const { return ts.days();         }     
00521     unsigned int day_fraction() const { return ts.day_fraction(); }
00522     int          sign()         const { return ts.sign();         }
00523     
00524     TimeStep GetTimeStep() const { return ts; };
00525     
00526     void     SetTimeStep(const TimeStep &ts_in) { ts = ts_in; }
00527     
00528     void SetDouble(const double d) { dts = d; }
00529     
00530   protected:
00531     TimeStep  ts;    
00532     double   dts;
00533   };
00534   
00535   UniverseTypeAwareTimeStep operator * (const int, const UniverseTypeAwareTimeStep &);
00536   UniverseTypeAwareTimeStep operator * (const UniverseTypeAwareTimeStep &, const int);
00537   
00538   UniverseTypeAwareTimeStep operator * (const double, const UniverseTypeAwareTimeStep &);
00539   UniverseTypeAwareTimeStep operator * (const UniverseTypeAwareTimeStep &, const double);
00540   
00541   // Angle
00542   
00543   class Angle {
00544   public:
00545     Angle() { radians = 0; }
00546     Angle(double x) { radians = x; }
00547     
00548   public:
00549     void   SetRad(double);
00550     void   GetRad(double&) const;
00551     double GetRad() const;
00552     //
00553     void   SetDPS(double  d, double  p, double  s);
00554     void   GetDPS(double &d, double &p, double &s) const;
00555     //
00556     void   SetHMS(double  h, double  m, double  s);
00557     void   GetHMS(double &h, double &m, double &s) const;
00558 
00559   private:
00560     double radians;
00561   };
00562   
00563   inline double sin(const Angle & alpha) {
00564     return std::sin(alpha.GetRad());
00565   }
00566   
00567   inline double cos(const Angle & alpha) {
00568     return std::cos(alpha.GetRad());
00569   }
00570   
00571   inline double tan(const Angle & alpha) {
00572     return std::tan(alpha.GetRad());
00573   }
00574   
00575   inline void sincos(const Angle & alpha, double & s, double & c) {
00576 #ifdef HAVE_SINCOS
00577     ::sincos(alpha.GetRad(),&s,&c); 
00578 #else // HAVE_SINCOS
00579     s = std::sin(alpha.GetRad());
00580     c = std::cos(alpha.GetRad());
00581 #endif // ! HAVE_SINCOS
00582   }
00583   
00584   // reference systems
00585   
00586   enum ReferenceSystem {
00587     EQUATORIAL=1,
00588     ECLIPTIC=2
00589   };
00590   
00591   inline void convert(ReferenceSystem &rs, const unsigned int i)  {
00592     switch(i) {
00593     case 1: rs = EQUATORIAL; break;
00594     case 2: rs = ECLIPTIC; break;
00595       //
00596     default:
00597       ORSA_ERROR("conversion problem: i = %i",i);
00598       break;
00599     }
00600   }
00601   
00602   Angle obleq(const Date &);
00603   // Angle delta_obleq(const Date&);
00604   // Angle delta_longitude(const Date&);
00605   Angle gmst(const Date &);
00606   
00607   // rotations  
00608   void EclipticToEquatorial(Vector&, const Date&);
00609   void EquatorialToEcliptic(Vector&, const Date&);
00610   
00611   // Angle gmst(const Date&);
00612   // void EquatorialToEcliptic(Vector&, const Date&);
00613   //
00614   // void EquatorialToDate_To_EquatorialJ2000(Vector&, const Date&);
00615   
00616   // J2000 versions
00617   Angle obleq_J2000();
00618   // Angle delta_obleq_J2000();
00619   // Angle delta_longitude_J2000();
00620   //
00621   void EclipticToEquatorial_J2000(Vector&);
00622   void EquatorialToEcliptic_J2000(Vector&);
00623   
00624 } // namespace orsa
00625 
00626 #include "orsa_file_jpl.h"
00627 
00628 namespace orsa {
00629   
00630   // another, improved orientation/rotation
00631   // formula, from the "Report of the IAU/IAG working group on
00632   // cartographic coordinates and rotational elements of the planets
00633   // and satellites: 2000", P.K. Seidelmann et al,
00634   // Celestial Mechanics and Dynamical Astronomy 82: 83-110 (2002).
00635   // data from Table I, pag. 86.
00636   void alpha_delta_meridian(const JPL_planets, const Date &,
00637                             Angle & alpha_zero, Angle & delta_zero, Angle & W);
00638   
00639 } // namespace orsa
00640 
00641 #endif // _ORSA_UNITS_H

Generated on Wed May 30 13:04:53 2007 for liborsa by  doxygen 1.5.2