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

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