orsa_units.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_units.h"
00026 
00027 #include "sdncal.h"
00028 
00029 #include "orsa_common.h"
00030 #include "orsa_universe.h"
00031 #include "orsa_file.h"
00032 
00033 #include <vector>
00034 #include <string>
00035 #include <iostream>
00036 
00037 #include <time.h>
00038 
00039 #include "support.h"
00040 
00041 using namespace std;
00042 
00043 namespace orsa {
00044   
00045   const double        G_MKS = 6.67259e-11;
00046   const double     MSUN_MKS = 1.9889e30;
00047   const double MJUPITER_MKS = 1.8989e27;
00048   const double   MEARTH_MKS = 5.9742e24;
00049   const double    MMOON_MKS = 7.3483e22;
00050   const double       AU_MKS = 1.49597870660e11;
00051   const double        c_MKS = 299792458.0;
00052   const double  R_EARTH_MKS = 6378140.0;
00053   const double   R_MOON_MKS = 1737400.0;
00054   
00055   double Units::GetG_MKS() const { return G_MKS; }
00056   
00057   Units::Units() {
00058     init_base();
00059     Time.Set(SECOND);
00060     Length.Set(M);
00061     Mass.Set(KG);
00062     TryToSetUnitsFromJPLFile();
00063     Recompute();
00064   };
00065   
00066   Units::Units(time_unit tu, length_unit lu, mass_unit mu) {
00067     init_base();
00068     Time.Set(tu);
00069     Length.Set(lu);
00070     Mass.Set(mu);
00071     TryToSetUnitsFromJPLFile();
00072     Recompute();
00073   }
00074   
00075   void Units::init_base() {
00076     G_base        =        G_MKS;
00077     MSun_base     =     MSUN_MKS;
00078     MJupiter_base = MJUPITER_MKS;
00079     MEarth_base   =   MEARTH_MKS;
00080     MMoon_base    =    MMOON_MKS;
00081     AU_base       =       AU_MKS;
00082     c_base        =        c_MKS;
00083     r_earth_base  =  R_EARTH_MKS;
00084     r_moon_base   =   R_MOON_MKS;
00085   }
00086   
00087   void Units::SetSystem(time_unit tu, length_unit lu, mass_unit mu) {
00088     Time.Set(tu);
00089     Length.Set(lu);
00090     Mass.Set(mu);
00091     Recompute();
00092   }
00093   
00094   void Units::TryToSetUnitsFromJPLFile() {
00095     
00096     /* 
00097        Config conf;
00098        OrsaConfigFile ocf(&conf);
00099        ocf.Read();
00100        ocf.Close();
00101     */
00102     if (config->paths[JPL_EPHEM_FILE]->GetValue() != "") {
00103       
00104       // cerr << "using JPL constants for units conversions..." << endl;
00105       
00106       JPLFile jf(config->paths[JPL_EPHEM_FILE]->GetValue().c_str());
00107       
00108       AU_base       = jf.GetAU_MKS();
00109       MSun_base     = jf.GetMSun_MKS();
00110       MJupiter_base = jf.GetMJupiter_MKS();
00111       MEarth_base   = jf.GetMEarth_MKS();
00112       MMoon_base    = jf.GetMMoon_MKS();
00113       c_base        = jf.GetC_MKS();
00114       r_earth_base  = jf.GetREarth_MKS();
00115       r_moon_base   = jf.GetRMoon_MKS();
00116       
00117     }
00118     
00119     Recompute();
00120   }
00121   
00122   double Units::GetTimeScale(const time_unit tu) const {
00123     if (tu == YEAR)   return (31557600.0);
00124     if (tu == DAY)    return (86400.0);
00125     if (tu == HOUR)   return (3600.0);
00126     if (tu == MINUTE) return (60.0);
00127     if (tu == SECOND) return (1.0);
00128     
00129     // here should throw an error...
00130     return 1.0;
00131   }
00132   
00133   string Units::label(const time_unit tu) const {
00134     if (tu == YEAR)   return "y";
00135     if (tu == DAY)    return "d";
00136     if (tu == HOUR)   return "h";
00137     if (tu == MINUTE) return "m";
00138     if (tu == SECOND) return "s";
00139     return "";
00140   }
00141   
00142   double Units::GetLengthScale(const length_unit lu) const {
00143     
00144     double ls=1.0;
00145     
00146     switch(lu) {
00147     case   MPARSEC: ls = (parsec_base*1.0e6); break;
00148     case   KPARSEC: ls = (parsec_base*1.0e3); break;
00149     case    PARSEC: ls = (parsec_base);       break;
00150     case        LY: ls = (c_base*31557600.0); break;
00151     case        AU: ls = (AU_base);           break;
00152     case EARTHMOON: ls = (3.844e8);           break;
00153     case    REARTH: ls = (r_earth_base);      break;
00154     case     RMOON: ls = (r_moon_base);       break;
00155     case        KM: ls = (1000.0);            break;
00156     case         M: ls = (1.0);               break;
00157     case        CM: ls = (0.01);              break;
00158     }
00159     
00160     return (ls);
00161   }
00162   
00163   string Units::label(const length_unit lu) const {
00164     if (lu == MPARSEC)   return "Mpc";
00165     if (lu == KPARSEC)   return "kpc";
00166     if (lu == PARSEC)    return "pc";
00167     if (lu == LY)        return "ly";
00168     if (lu == AU)        return "AU";
00169     if (lu == EARTHMOON) return "LD"; // earth-moon mean distance
00170     if (lu == REARTH)    return "ER";
00171     if (lu == RMOON)     return "MR";
00172     if (lu == KM)        return "km";
00173     if (lu ==  M)        return "m";
00174     if (lu == CM)        return "cm";
00175     return "";
00176   }
00177   
00178   double Units::GetMassScale(const mass_unit mu) const {
00179     
00180     double ms=1.0;
00181     
00182     switch(mu) {
00183     case     MSUN: ms = (MSun_base);     break;
00184     case MJUPITER: ms = (MJupiter_base); break;
00185     case   MEARTH: ms = (MEarth_base);   break;
00186     case    MMOON: ms = (MMoon_base);    break;
00187     case       KG: ms = (1.0);           break;
00188     case     GRAM: ms = (0.001);         break;
00189     }
00190     
00191     return (ms);
00192   }
00193       
00194   string Units::label(const mass_unit mu) const {
00195     if (mu == MSUN)     return "Sun mass";
00196     if (mu == MJUPITER) return "Jupiter mass";
00197     if (mu == MEARTH)   return "Earth mass";
00198     if (mu == MMOON)    return "Moon mass";
00199     if (mu == KG)       return "kg";
00200     if (mu == GRAM)     return "g";
00201     return "";
00202   }
00203   
00204   double Units::GetTimeScale()   const { return Units::GetTimeScale(     Time.GetBaseUnit()); } 
00205   double Units::GetLengthScale() const { return Units::GetLengthScale( Length.GetBaseUnit()); } 
00206   double Units::GetMassScale()   const { return Units::GetMassScale(     Mass.GetBaseUnit()); } 
00207   
00208   void Units::Recompute() {
00209     G = G_base*secure_pow(GetTimeScale(),2.0)*secure_pow(GetLengthScale(),-3.0)*secure_pow(GetMassScale(),1.0); 
00210     MSun = MSun_base/GetMassScale();
00211     c = c_base*GetTimeScale()/GetLengthScale();
00212     parsec_base = AU_base/(2*sin((pi/180)/3600.0/2));
00213   }
00214   
00215   // Date related tables
00216   
00217   // const TimeScale default_Date_timescale = TDT; // ET
00218   TimeScale default_Date_timescale = TT; // ET, updated by the Universe costructor
00219   
00220   struct TAI_minus_UTC_element {
00221     int day,month,year;
00222     int TAI_minus_UTC;
00223   };
00224   
00225   inline bool operator == (const TAI_minus_UTC_element &x, const TAI_minus_UTC_element &y) {
00226     if (x.day   != y.day)   return false;
00227     if (x.month != y.month) return false;
00228     if (x.year  != y.year)  return false;
00229     if (x.TAI_minus_UTC != y.TAI_minus_UTC) return false;
00230     return true;
00231   }
00232   
00233   inline bool operator != (const TAI_minus_UTC_element &x, const TAI_minus_UTC_element &y) {
00234     return (!(x==y));
00235   }
00236   
00237   /*
00238    *  The TAI_minus_UTC_table should be updated periodically, 
00239    * for instance by consulting:
00240    * 
00241    * http://hpiers.obspm.fr/webiers/general/communic/publi/PUBLI.html
00242    * 
00243    * The present version is updated through Bulletin C17 (Jan 28, 1999)
00244    * Day, month, year, TAI-UTC (s)
00245    */
00246   
00247   const TAI_minus_UTC_element TAI_minus_UTC_table_final_element = {0,0,0,0};
00248   const TAI_minus_UTC_element TAI_minus_UTC_table[] = {
00249     {1,1,1972,10},
00250     {1,7,1972,11},
00251     {1,1,1973,12},
00252     {1,1,1974,13},
00253     {1,1,1975,14},
00254     {1,1,1976,15},
00255     {1,1,1977,16},
00256     {1,1,1978,17},
00257     {1,1,1979,18},
00258     {1,1,1980,19},
00259     {1,7,1981,20},
00260     {1,7,1982,21},
00261     {1,7,1983,22},
00262     {1,7,1985,23},
00263     {1,1,1988,24},
00264     {1,1,1990,25},
00265     {1,1,1991,26},
00266     {1,7,1992,27},
00267     {1,7,1993,28},
00268     {1,7,1994,29},
00269     {1,1,1996,30},
00270     {1,7,1997,31},
00271     {1,1,1999,32}, // latest leap second, updated Jan 2003
00272     {0,0,0,0} // TAI_minus_UTC_table_final_element
00273   };
00274   
00275   struct ET_minus_UT_element {
00276     int day,month,year;
00277     double ET_minus_UT;
00278   };
00279   
00280   inline bool operator == (const ET_minus_UT_element &x, const ET_minus_UT_element &y) {
00281     if (x.day   != y.day)   return false;
00282     if (x.month != y.month) return false;
00283     if (x.year  != y.year)  return false;
00284     if (x.ET_minus_UT != y.ET_minus_UT) return false;
00285     return true;
00286   }
00287   
00288   inline bool operator != (const ET_minus_UT_element &x, const ET_minus_UT_element &y) {
00289     return (!(x==y));
00290   }
00291   
00292   /* 
00293    * Values ET - UT at 0h UT of the date (specified as day, month, year)
00294    * from The Astronomical Almanac 1999, pages K8-9
00295    */
00296   
00297   const ET_minus_UT_element ET_minus_UT_table_final_element = {0,0,0,0};
00298   const ET_minus_UT_element ET_minus_UT_table[] = {
00299     {1,1,1800,13.7},
00300     {1,1,1801,13.4},
00301     {1,1,1802,13.1},
00302     {1,1,1803,12.9},
00303     {1,1,1804,12.7},
00304     {1,1,1805,12.6},
00305     {1,1,1806,12.5},
00306     {1,1,1816,12.5},
00307     {1,1,1817,12.4},
00308     {1,1,1818,12.3},
00309     {1,1,1819,12.2},
00310     {1,1,1820,12.0},
00311     {1,1,1821,11.7},
00312     {1,1,1822,11.4},
00313     {1,1,1823,11.1},
00314     {1,1,1824,10.6},
00315     {1,1,1825,10.2},
00316     {1,1,1826, 9.6},
00317     {1,1,1827, 9.1},
00318     {1,1,1828, 8.6},
00319     {1,1,1829, 8.0},
00320     {1,1,1830, 7.5},
00321     {1,1,1831, 7.0},
00322     {1,1,1832, 6.6},
00323     {1,1,1833, 6.3},
00324     {1,1,1834, 6.0},
00325     {1,1,1835, 5.8},
00326     {1,1,1836, 5.7},
00327     {1,1,1837, 5.6},
00328     {1,1,1838, 5.6},
00329     {1,1,1839, 5.6},
00330     {1,1,1840, 5.7},
00331     {1,1,1841, 5.8},
00332     {1,1,1842, 5.9},
00333     {1,1,1843, 6.1},
00334     {1,1,1844, 6.2},
00335     {1,1,1845, 6.3},
00336     {1,1,1846, 6.5},
00337     {1,1,1847, 6.6},
00338     {1,1,1848, 6.8},
00339     {1,1,1849, 6.9},
00340     {1,1,1850, 7.1},
00341     {1,1,1851, 7.2},
00342     {1,1,1852, 7.3},
00343     {1,1,1853, 7.4},
00344     {1,1,1854, 7.5},
00345     {1,1,1855, 7.6},
00346     {1,1,1856, 7.7},
00347     {1,1,1857, 7.7},
00348     {1,1,1858, 7.8},
00349     {1,1,1859, 7.8},
00350     {1,1,1860, 7.88},
00351     {1,1,1861, 7.82},
00352     {1,1,1862, 7.54},
00353     {1,1,1863, 6.97},
00354     {1,1,1864, 6.40},
00355     {1,1,1865, 6.02},
00356     {1,1,1866, 5.41},
00357     {1,1,1867, 4.10},
00358     {1,1,1868, 2.92},
00359     {1,1,1869, 1.82},
00360     {1,1,1870, 1.61},
00361     {1,1,1871, 0.10},
00362     {1,1,1872,-1.02},
00363     {1,1,1873,-1.28},
00364     {1,1,1874,-2.69},
00365     {1,1,1875,-3.24},
00366     {1,1,1876,-3.64},
00367     {1,1,1877,-4.54},
00368     {1,1,1878,-4.71},
00369     {1,1,1879,-5.11},
00370     {1,1,1880,-5.40},
00371     {1,1,1881,-5.42},
00372     {1,1,1882,-5.20},
00373     {1,1,1883,-5.46},
00374     {1,1,1884,-5.46},
00375     {1,1,1885,-5.79},
00376     {1,1,1886,-5.63},
00377     {1,1,1887,-5.64},
00378     {1,1,1888,-5.80},
00379     {1,1,1889,-5.66},
00380     {1,1,1890,-5.87},
00381     {1,1,1891,-6.01},
00382     {1,1,1892,-6.19},
00383     {1,1,1893,-6.64},
00384     {1,1,1894,-6.44},
00385     {1,1,1895,-6.47},
00386     {1,1,1896,-6.09},
00387     {1,1,1897,-5.76},
00388     {1,1,1898,-4.66},
00389     {1,1,1899,-3.74},
00390     {1,1,1900,-2.72},
00391     {1,1,1901,-1.54},
00392     {1,1,1902,-0.02},
00393     {1,1,1903, 1.24},
00394     {1,1,1904, 2.64},
00395     {1,1,1905, 3.86},
00396     {1,1,1906, 5.37},
00397     {1,1,1907, 6.14},
00398     {1,1,1908, 7.75},
00399     {1,1,1909, 9.13},
00400     {1,1,1910,10.46},
00401     {1,1,1911,11.53},
00402     {1,1,1912,13.36},
00403     {1,1,1913,14.65},
00404     {1,1,1914,16.01},
00405     {1,1,1915,17.20},
00406     {1,1,1916,18.24},
00407     {1,1,1917,19.06},
00408     {1,1,1918,20.25},
00409     {1,1,1919,20.95},
00410     {1,1,1920,21.16},
00411     {1,1,1921,22.25},
00412     {1,1,1922,22.41},
00413     {1,1,1923,23.03},
00414     {1,1,1924,23.49},
00415     {1,1,1925,23.62},
00416     {1,1,1926,23.86},
00417     {1,1,1927,24.49},
00418     {1,1,1928,24.34},
00419     {1,1,1929,24.08},
00420     {1,1,1930,24.02},
00421     {1,1,1931,24.00},
00422     {1,1,1932,23.87},
00423     {1,1,1933,23.95},
00424     {1,1,1934,23.86},
00425     {1,1,1935,23.93},
00426     {1,1,1936,23.73},
00427     {1,1,1937,23.92},
00428     {1,1,1938,23.96},
00429     {1,1,1939,24.02},
00430     {1,1,1940,24.33},
00431     {1,1,1941,24.83},
00432     {1,1,1942,25.30},
00433     {1,1,1943,25.70},
00434     {1,1,1944,26.24},
00435     {1,1,1945,26.77},
00436     {1,1,1946,27.28},
00437     {1,1,1947,27.78},
00438     {1,1,1948,28.25},
00439     {1,1,1949,28.71},
00440     {1,1,1950,29.15},
00441     {1,1,1951,29.57},
00442     {1,1,1952,29.97},
00443     {1,1,1953,30.36},
00444     {1,1,1954,30.72},
00445     {1,1,1955,31.07},
00446     {1,1,1956,31.35},
00447     {1,1,1957,31.68},
00448     {1,1,1958,32.18},
00449     {1,1,1959,32.68},
00450     {1,1,1960,33.15},
00451     {1,1,1961,33.59},
00452     {1,1,1962,34.00},
00453     {1,1,1963,34.47},
00454     {1,1,1964,35.03},
00455     {1,1,1965,35.73},
00456     {1,1,1966,36.54},
00457     {1,1,1967,37.43},
00458     {1,1,1968,38.29},
00459     {1,1,1969,39.20},
00460     {1,1,1970,40.18},
00461     {1,1,1971,41.17},
00462     {1,1,1972,42.23},
00463     {1,1,1973,43.37},
00464     {1,1,1974,44.49},
00465     {1,1,1975,45.48},
00466     {1,1,1976,46.46},
00467     {1,1,1977,47.52},
00468     {1,1,1978,48.53},
00469     {1,1,1979,49.59},
00470     {1,1,1980,50.54},
00471     {1,1,1981,51.38},
00472     {1,1,1982,52.17},
00473     {1,1,1983,52.96},
00474     {1,1,1984,53.79},
00475     {1,1,1985,54.34},
00476     {1,1,1986,54.87},
00477     {1,1,1987,55.32},
00478     {1,1,1988,55.82},
00479     {1,1,1989,56.30},
00480     {1,1,1990,56.86},
00481     {1,1,1991,57.57},
00482     {1,1,1992,58.31},
00483     {1,1,1993,59.12},
00484     {1,1,1994,59.98},
00485     {1,1,1995,60.78},
00486     {1,1,1996,61.63},
00487     {1,1,1997,62.29},
00488     {1,1,1998,62.97},
00489     {1,1,1999,63.46},
00490     {1,1,2000,63.83},
00491     {1,1,2001,64.09},
00492     {1,1,2002,64.30},
00493     {1,1,2003,64.47},
00494     {0,0,0,0} // ET_minus_UT_table_final_element
00495   };
00496   
00497   // Date
00498   
00499   double Date::delta_seconds(int y, int m, int d, const TimeScale from, const TimeScale to) {
00500     
00501     double delta_seconds=0.0;
00502     
00503     if (to==from) return delta_seconds;
00504     
00505     int j = 0;
00506     
00507     // convert from 'from' to TDT
00508     switch (from) {
00509       
00510     case TDT: 
00511       delta_seconds=0.0;
00512       break;
00513       
00514     case TAI:
00515       delta_seconds=32.184;
00516       break;
00517       
00518     case UTC:
00519       delta_seconds = 32.184;
00520       if (y>=TAI_minus_UTC_table[0].year) {
00521         j=0;
00522         TAI_minus_UTC_element candidate = TAI_minus_UTC_table[0];
00523         while (TAI_minus_UTC_table[j] != TAI_minus_UTC_table_final_element) {
00524           
00525           if (TAI_minus_UTC_table[j].year   <=y) {
00526             if (TAI_minus_UTC_table[j].month<=m) {
00527               if (TAI_minus_UTC_table[j].day<=d) {
00528                 if ( ( (TAI_minus_UTC_table[j].year >  candidate.year) ) ||
00529                      ( (TAI_minus_UTC_table[j].year == candidate.year) && (TAI_minus_UTC_table[j].month >  candidate.month) ) ||
00530                      ( (TAI_minus_UTC_table[j].year == candidate.year) && (TAI_minus_UTC_table[j].month == candidate.month) && (TAI_minus_UTC_table[j].day > candidate.day) ) ) { 
00531                   candidate = TAI_minus_UTC_table[j];
00532                 }               
00533               }
00534             }
00535           }
00536           
00537           j++;
00538         }
00539         
00540         delta_seconds += candidate.TAI_minus_UTC;
00541       }
00542       break;
00543       
00544     case UT1:
00545       delta_seconds = 0.0;
00546       if (y>=ET_minus_UT_table[0].year) {
00547         j=0;
00548         ET_minus_UT_element candidate=ET_minus_UT_table[0];
00549         while (ET_minus_UT_table[j] != ET_minus_UT_table_final_element) {
00550           
00551           if (ET_minus_UT_table[j].year   <=y) {
00552             if (ET_minus_UT_table[j].month<=m) {
00553               if (ET_minus_UT_table[j].day<=d) {
00554                 
00555                 if ( ( (ET_minus_UT_table[j].year >  candidate.year) ) ||
00556                      ( (ET_minus_UT_table[j].year == candidate.year) && (ET_minus_UT_table[j].month >  candidate.month) ) ||
00557                      ( (ET_minus_UT_table[j].year == candidate.year) && (ET_minus_UT_table[j].month == candidate.month) && (ET_minus_UT_table[j].day > candidate.day) ) ) {
00558                   candidate = ET_minus_UT_table[j];
00559                 }
00560               }
00561             }
00562           }
00563           
00564           j++;
00565         }
00566         delta_seconds += candidate.ET_minus_UT;
00567       }
00568       break;
00569       
00570     case GPS:
00571       delta_seconds = 32.184 + 19.0;
00572       break;
00573     }
00574     
00575     // convert from TDT to 'to'
00576     switch (to) {
00577       
00578     case TDT:
00579       break;
00580       
00581     case TAI: 
00582       delta_seconds -= 32.184; 
00583       break;
00584       
00585     case UTC: 
00586       delta_seconds -= 32.184;
00587       if (y>=TAI_minus_UTC_table[0].year) {
00588         j=0;
00589         TAI_minus_UTC_element candidate = TAI_minus_UTC_table[0];
00590         while (TAI_minus_UTC_table[j] != TAI_minus_UTC_table_final_element) {
00591           
00592           if (TAI_minus_UTC_table[j].year   <=y) {
00593             if (TAI_minus_UTC_table[j].month<=m) {
00594               if (TAI_minus_UTC_table[j].day<=d) {
00595                 if ( ( (TAI_minus_UTC_table[j].year >  candidate.year) ) ||
00596                      ( (TAI_minus_UTC_table[j].year == candidate.year) && (TAI_minus_UTC_table[j].month >  candidate.month) ) ||
00597                      ( (TAI_minus_UTC_table[j].year == candidate.year) && (TAI_minus_UTC_table[j].month == candidate.month) && (TAI_minus_UTC_table[j].day > candidate.day) ) ) { 
00598                   candidate = TAI_minus_UTC_table[j];
00599                 }
00600               }
00601             }
00602           }
00603           
00604           j++;
00605         }
00606         
00607         delta_seconds -= candidate.TAI_minus_UTC;
00608       }  
00609       break;
00610       
00611     case UT1:
00612       if (y>=ET_minus_UT_table[0].year) {
00613         j=0;
00614         ET_minus_UT_element candidate=ET_minus_UT_table[0];
00615         while (ET_minus_UT_table[j] != ET_minus_UT_table_final_element) {
00616           
00617           if (ET_minus_UT_table[j].year   <=y) {
00618             if (ET_minus_UT_table[j].month<=m) {
00619               if (ET_minus_UT_table[j].day<=d) {
00620                 if ( ( (ET_minus_UT_table[j].year >  candidate.year) ) ||
00621                      ( (ET_minus_UT_table[j].year == candidate.year) && (ET_minus_UT_table[j].month >  candidate.month) ) ||
00622                      ( (ET_minus_UT_table[j].year == candidate.year) && (ET_minus_UT_table[j].month == candidate.month) && (ET_minus_UT_table[j].day > candidate.day) ) ) {
00623                   candidate = ET_minus_UT_table[j];
00624                 }
00625               }
00626             }
00627           }
00628           
00629           j++;
00630         }
00631         
00632         delta_seconds -= candidate.ET_minus_UT;
00633       }
00634       break;
00635       
00636     case GPS: 
00637       delta_seconds -= (32.184 + 19.0);
00638       break;
00639     }
00640     
00641     // cerr << "delta_seconds: from " << TimeScaleLabel(from) << " to " << TimeScaleLabel(to) << " = " << delta_seconds << endl;
00642     
00643     return (delta_seconds);
00644   }
00645   
00646   Date::Date() : sdn(0), df(0) { }
00647   
00648   Date::Date(const UniverseTypeAwareTime & t) : sdn(t.GetDate().sdn), df(t.GetDate().df) { }
00649   
00650   Date::Date(const Date & d) : sdn(d.sdn), df(d.df) { }
00651   
00652   //! fractionary day, relative to 00:00  
00653   void Date::SetGregor(int  y, int  m, double d, TimeScale ts) {
00654     
00655     int id = (int)(floor(d));
00656     d -= id;
00657     //
00658     d *= 24;
00659     int H = (int)(floor(d));
00660     d -= H;
00661     //
00662     d *= 60;
00663     int M = (int)(floor(d));
00664     d -= M;
00665     //
00666     d *= 60;
00667     int S = (int)(floor(d));
00668     d -= S;
00669     //
00670     d *= 1000;
00671     int ms = (int)(floor(d));
00672     
00673     SetGregor(y,m,id,H,M,S,ms,ts);
00674   }
00675   
00676   void Date::SetGregor(int y, int m, int d, int H, int M, int S, int ms, TimeScale ts) {
00677     
00678     ms += (int)(1000.0 * delta_seconds(y,m,d,ts));
00679     
00680     // checks...
00681     while (ms >= 1000) {
00682       S  += 1;
00683       ms -= 1000;
00684     }
00685     while (S >= 60) {
00686       M += 1;
00687       S -= 60;
00688     }
00689     while (M >= 60) {
00690       H += 1;
00691       M -= 60;
00692     }
00693     while (H >= 24) {
00694       d += 1;
00695       H -= 24;
00696     }
00697     
00698     // more...
00699     while (ms < 0) {
00700       S  -= 1;
00701       ms += 60;
00702     }
00703     while (S < 0) {
00704       M -= 1;
00705       S += 60;
00706     }
00707     while (M < 0) {
00708       H -= 1;
00709       M += 60;
00710     }
00711     while (H < 0) {
00712       d -= 1;
00713       H += 24;
00714     }
00715     
00716     sdn = GregorianToSdn(y,m,d);
00717     df  = ( ( (H * 60 + M) * 60 + S) * 1000 + ms) * (TimeStep::max_day_fraction() / 86400000);
00718   }
00719   
00720   /* 
00721      void Date::GetGregor(int &y, int &m, int &d, TimeScale ts) const {
00722      SdnToGregorian((long int)floor(_ref_time+_delta_time), &y, &m, &d);
00723      const double tmp1 = _ref_time   - floor(_ref_time);
00724      const double tmp2 = _delta_time - floor(_delta_time);
00725      const double local_frac = (tmp1+tmp2) - floor(tmp1+tmp2) - (1.0/86400.0)*delta_seconds(y,m,d,ts); 
00726      if (local_frac < 0.0) {
00727      SdnToGregorian((long int)floor(_ref_time+_delta_time)-1, &y, &m, &d);
00728      } else if (local_frac > 1.0) { 
00729      SdnToGregorian((long int)floor(_ref_time+_delta_time)+1, &y, &m, &d);
00730      }
00731      }
00732   */
00733   
00734   void Date::GetGregor(int &y, int &m, int &d, TimeScale ts) const {
00735     SdnToGregorian(sdn, &y, &m, &d);
00736     const double ds = delta_seconds(y,m,d,ts);
00737     const int delta_df = - (int)(ds*(TimeStep::max_day_fraction() / 86400)); // negative!
00738     if (delta_df < 0) {
00739       if (abs(delta_df) > df) {
00740         SdnToGregorian(sdn-1, &y, &m, &d);
00741       }
00742     } else {
00743       const unsigned int local_df = df + delta_df;
00744       if (local_df >= TimeStep::max_day_fraction()) {
00745         SdnToGregorian(sdn+1, &y, &m, &d);
00746       }
00747     }
00748   }
00749   
00750   double Date::GetDayFraction(TimeScale ts) const {
00751     return (GetDayFraction_unsigned_int(ts)*1.0)/(TimeStep::max_day_fraction());
00752   }
00753   
00754   unsigned int Date::GetDayFraction_unsigned_int(TimeScale ts) const {
00755     int y,m,d;
00756     SdnToGregorian(sdn, &y, &m, &d);
00757     const double ds = delta_seconds(y,m,d,ts);
00758     const int delta_df = - (int)(ds*(TimeStep::max_day_fraction() / 86400)); // negative!
00759     unsigned int local_df = 0;
00760     if (delta_df < 0) {
00761       if (abs(delta_df) > df) {
00762         local_df = TimeStep::max_day_fraction() + df - abs(delta_df);
00763       } else {
00764         local_df = df + delta_df;
00765       }
00766     } else {
00767       local_df = df + delta_df;
00768     }
00769     local_df %= TimeStep::max_day_fraction();
00770     return (local_df);
00771   }
00772   
00773   double Date::Time() const {
00774     return (FromUnits(GetJulian(),DAY));
00775   }
00776   
00777   double Date::GetTime() const {
00778     return (FromUnits(GetJulian(),DAY));
00779   }
00780   
00781   void Date::SetTime(const double t) {
00782     SetJulian(FromUnits(t,DAY,-1));
00783   }
00784   
00785   Date & Date::operator += (const UniverseTypeAwareTimeStep & ts) {
00786     
00787     sdn += ts.sign() * ts.days();
00788     
00789     if (ts.sign() == -1) {
00790       if (ts.day_fraction() > df) {
00791         --sdn;
00792         df = TimeStep::max_day_fraction() + df - ts.day_fraction();
00793       } else {
00794         df -= ts.day_fraction();
00795       }
00796     } else {
00797       df += ts.day_fraction();
00798     }
00799     
00800     while (df >= TimeStep::max_day_fraction()) {
00801       ++sdn;
00802       df -= TimeStep::max_day_fraction();
00803     }
00804     
00805     return * this;
00806   }
00807   
00808   Date & Date::operator -= (const UniverseTypeAwareTimeStep & ts) {
00809     * this += - ts;
00810     return * this;
00811   }
00812   
00813   bool Date::operator == (const Date & d) const {
00814     if (sdn != d.sdn) return false;
00815     if (df != d.df) return false;
00816     return true;
00817   }
00818   
00819   bool Date::operator != (const Date & d) const {
00820     return (!((*this)==d));
00821   }
00822   
00823   bool Date::operator < (const Date & d) const {
00824     return (GetTimeStep() < d.GetTimeStep());
00825   }
00826   
00827   bool Date::operator > (const Date & d) const {
00828     return (GetTimeStep() > d.GetTimeStep());
00829   }
00830   
00831   void Date::SetJulian(double jd, TimeScale ts) {
00832     sdn = (unsigned int)(floor(jd));
00833     double frac = jd - floor(jd) + 0.5;
00834     if (frac >= 1.0) {
00835       ++sdn;
00836       frac = fmod(frac,1.0);
00837     }
00838     //
00839     {
00840       int y,m,d;
00841       SdnToGregorian(sdn, &y, &m, &d);
00842       jd += (1.0/86400.0)*delta_seconds(y,m,d,ts);
00843     }
00844     sdn = (unsigned int)(floor(jd));
00845     frac = jd - floor(jd) + 0.5;
00846     if (frac >= 1.0) {
00847       ++sdn;
00848       frac = fmod(frac,1.0);
00849     }
00850     
00851     df = (unsigned int)rint(frac * TimeStep::max_day_fraction());
00852   }
00853   
00854   void Date::GetJulian(double & jd, TimeScale ts) const {
00855     int y,m,d;
00856     SdnToGregorian(sdn, &y, &m, &d);
00857     jd = sdn + ((df*1.0)/ TimeStep::max_day_fraction()) - 0.5 - (1.0/86400.0)*delta_seconds(y,m,d,ts);
00858   }
00859   
00860   double Date::GetJulian(TimeScale ts) const {
00861     double jd;
00862     GetJulian(jd,ts);
00863     return jd;
00864   }
00865   
00866   void Date::SetNow() {
00867     const time_t tt_now = time(0);
00868     struct tm *tm_struct = gmtime(&tt_now);
00869     // this is accurate down to a second, microseconds are not used
00870     // if they were, the term is ...tm_sec+(tv.tv_usec*1.0e-6))...
00871     // SetGregor(1900+tm_struct->tm_year,1+tm_struct->tm_mon,tm_struct->tm_mday+((tm_struct->tm_sec)/86400.0+tm_struct->tm_min/1440.0+tm_struct->tm_hour/24.0),UTC);
00872     SetGregor(1900+tm_struct->tm_year,
00873               1+tm_struct->tm_mon,
00874               tm_struct->tm_mday,
00875               tm_struct->tm_hour,
00876               tm_struct->tm_min,
00877               tm_struct->tm_sec,
00878               0, // ms
00879               UTC);
00880   }
00881   
00882   void Date::SetToday() {
00883     SetNow();
00884     df = 0; // reset day fraction
00885   }
00886   
00887   void Date::SetJ2000() {
00888     SetJulian(2451545.0,TT); /* IAU definition: [2000 Jan 1d 12h TDT]
00889                                 http://en.wikipedia.org/wiki/Terrestrial_Time 
00890                                 http://en.wikipedia.org/wiki/Month 
00891                                 http://aa.usno.navy.mil/software/novas/novas_c/novasc_info.html
00892                                 http://nanotitan.com/software/api/suite/diamond/nT/quantity/constant/TIME_INSTANT.html#J2000
00893                              */
00894   }
00895   
00896   string TimeScaleLabel(TimeScale ts) {
00897     if (ts==UTC) return "UTC";
00898     if (ts==UT)  return "UT";
00899     if (ts==UT1) return "UT1";
00900     if (ts==TAI) return "TAI";
00901     if (ts==TDT) return "TDT";
00902     if (ts==ET)  return "ET";
00903     if (ts==GPS) return "GPS";     
00904     return "";
00905   }
00906   
00907   // TimeStep
00908   
00909   TimeStep::TimeStep() : _days(0), _day_fraction(0), _sign(+1) { 
00910     internal_check();
00911   }
00912   
00913   TimeStep::TimeStep(const unsigned int days, const unsigned int day_fraction, const int sign) : _days(days), _day_fraction(day_fraction), _sign(sign) {
00914     if (_sign == 0) {
00915       ORSA_ERROR("Hmmm, sign equal to zero...");
00916     }   else {
00917       _sign = _sign / abs(_sign);
00918     }
00919     internal_check();
00920   }
00921   
00922   TimeStep::TimeStep(const double t) {
00923     if (t < 0) {
00924       _sign = -1;
00925     } else {
00926       _sign = +1;
00927     }   
00928     
00929     const double t_in_days = FromUnits(fabs(t),DAY,-1);
00930     _days         = (unsigned int)(floor(t_in_days));
00931     _day_fraction = (unsigned int)rint((t_in_days - _days)*max_day_fraction());
00932     
00933     internal_check();
00934   }
00935   
00936   TimeStep::TimeStep(const TimeStep & ts) : _days(ts._days), _day_fraction(ts._day_fraction), _sign(ts._sign) {
00937     internal_check();
00938   }
00939   
00940   void TimeStep::internal_check() {
00941     if (IsZero()) _sign = +1;
00942   }
00943   
00944   void TimeStep::AddDays(const unsigned int d, const int sign) {
00945     if (sign == _sign) {
00946       _days += d;
00947     } else {
00948       if (d > _days) {
00949         _sign = -_sign;
00950         _days = d - _days - 1;
00951         _day_fraction = max_day_fraction() - _day_fraction;
00952         if (_day_fraction >= max_day_fraction()) {
00953           ++_days;
00954           _day_fraction -= max_day_fraction();
00955         }
00956       } else {
00957         _days -= d;
00958       }
00959     }
00960     internal_check();
00961   }
00962   
00963   void TimeStep::AddDayFractions(const unsigned int df, const int sign) {
00964     if (sign == _sign) {
00965       _day_fraction += df;
00966       if (_day_fraction >= max_day_fraction()) {
00967         ++_days;
00968         _day_fraction -= max_day_fraction();
00969       }
00970     } else {
00971       if (_day_fraction >= df) {
00972         _day_fraction -= df;
00973       } else {
00974         if (_days > 0) {
00975           --_days;
00976           _day_fraction += max_day_fraction();
00977           _day_fraction -= df;
00978         } else {
00979           _sign = -_sign;
00980           _day_fraction = df - _day_fraction;
00981         }
00982       }
00983     }
00984     internal_check();
00985   }
00986   
00987   /* 
00988      TimeStep & TimeStep::operator += (const TimeStep & t) { 
00989      
00990      if (_sign != t._sign) {
00991      if ( (t._days > _days) || 
00992      ( (t._days == _days) && 
00993      (_day_fraction >= t._day_fraction) ) ) {
00994      _sign = -_sign;
00995      _days = t._days - _days;
00996      if (_day_fraction >= t._day_fraction) {
00997      _day_fraction -= t._day_fraction;
00998      } else {
00999      --_days;
01000      _day_fraction += max_day_fraction();
01001      _day_fraction -= t._day_fraction;
01002      }
01003      } else {
01004      _days  -= t._days;
01005      
01006      // _day_fraction -= t._day_fraction;
01007      if (_day_fraction >= t._day_fraction) {
01008      _day_fraction -= t._day_fraction;
01009      } else {
01010      --_days;
01011      _day_fraction += max_day_fraction();
01012      _day_fraction -= t._day_fraction;
01013      }
01014      
01015      }
01016      } else {
01017      _days         += t._days;
01018      _day_fraction += t._day_fraction;
01019      }
01020      
01021      while (_day_fraction >= max_day_fraction()) {
01022      ++_days;
01023      _day_fraction -= max_day_fraction();
01024      }
01025      
01026      return * this;
01027      }
01028   */
01029   
01030   /* 
01031      TimeStep & TimeStep::operator += (const TimeStep & ts) { 
01032      if (_sign == ts._sign) {
01033      _days         += ts._days;
01034      _day_fraction += ts._day_fraction;
01035      while (_day_fraction >= max_day_fraction()) {
01036      ++_days;
01037      _day_fraction -= max_day_fraction();
01038      }
01039      } else {
01040      * this -= - ts;
01041      }
01042      return * this;
01043      }
01044   */
01045   
01046   TimeStep & TimeStep::operator += (const TimeStep & ts) { 
01047     AddDays(ts._days,ts._sign);
01048     AddDayFractions(ts._day_fraction,ts._sign);
01049     return * this;
01050   }
01051   
01052   TimeStep & TimeStep::operator -= (const TimeStep & ts) { 
01053     AddDays(ts._days,-ts._sign);
01054     AddDayFractions(ts._day_fraction,-ts._sign);
01055     return * this;
01056   }
01057   
01058   /* 
01059      TimeStep & TimeStep::operator -= (const TimeStep & ts) { 
01060      if (_sign == ts._sign) {
01061      if (ts._days > _days) {
01062      _sign = -_sign;
01063      _days = ts._days - _days;
01064      if (_day_fraction >= ts._day_fraction) {
01065      _day_fraction -= ts._day_fraction;
01066      } else {
01067      --_days;
01068      _day_fraction += max_day_fraction();
01069      _day_fraction -= ts._day_fraction;
01070      }
01071      } else if ( (ts._days == _days) && 
01072      (ts._day_fraction > _day_fraction) ) {
01073      _sign = -_sign;
01074      _days = 0;
01075      _day_fraction = ts._day_fraction - _day_fraction;
01076      } else {
01077      _days -= ts._days;
01078      // _day_fraction -= ts._day_fraction;
01079      if (_day_fraction >= ts._day_fraction) {
01080      _day_fraction -= ts._day_fraction;
01081      } else {
01082      --_days;
01083      _day_fraction += max_day_fraction();
01084      _day_fraction -= ts._day_fraction;
01085      }
01086      }
01087      } else {
01088      * this += - ts;
01089      }      
01090      return * this;
01091      }
01092   */
01093   
01094   TimeStep TimeStep::operator + (const TimeStep & ts) const {
01095     TimeStep _ts(*this);
01096     _ts += ts;
01097     return _ts;
01098   }
01099   
01100   TimeStep TimeStep::operator - (const TimeStep & ts) const {
01101     TimeStep _ts(*this);
01102     _ts -= ts;
01103     return _ts;
01104   }
01105   
01106   /* 
01107      TimeStep & TimeStep::operator *= (const int p) {
01108      
01109      if (p == 0) {
01110      _days = _day_fraction = 0;
01111      _sign = +1;
01112      return * this;
01113      }
01114      
01115      if (p < 0) _sign = -_sign;
01116      
01117      _days *= abs(p);  
01118      
01119      const unsigned int original_day_fraction = _day_fraction;
01120      unsigned int count = abs(p)-1;
01121      while (count != 0) {
01122      _day_fraction += original_day_fraction;
01123      while (_day_fraction >= max_day_fraction()) {
01124      ++_days;
01125      _day_fraction -= max_day_fraction();
01126      }
01127      --count;
01128      }
01129      
01130      return * this;
01131      }
01132   */
01133   
01134   TimeStep & TimeStep::operator *= (const int p) {
01135     const double t = GetDouble()*p;
01136     // same code in constructor from double
01137     if (t < 0) {
01138       _sign = -1;
01139     } else {
01140       _sign = +1;
01141     }   
01142     //    
01143     const double t_in_days = FromUnits(fabs(t),DAY,-1);
01144     _days         = (unsigned int)(floor(t_in_days));
01145     _day_fraction = (unsigned int)rint((t_in_days - _days)*max_day_fraction());
01146     //
01147     internal_check();
01148     //
01149     return * this;
01150   }     
01151   
01152   TimeStep & TimeStep::operator *= (const double x) {
01153     const double t = GetDouble()*x;
01154     // same code in constructor from double
01155     if (t < 0) {
01156       _sign = -1;
01157     } else {
01158       _sign = +1;
01159     }   
01160     //    
01161     const double t_in_days = FromUnits(fabs(t),DAY,-1);
01162     _days         = (unsigned int)(floor(t_in_days));
01163     _day_fraction = (unsigned int)rint((t_in_days - _days)*max_day_fraction());
01164     //
01165     internal_check();
01166     //
01167     return * this;
01168   }
01169   
01170   /* 
01171      TimeStep & TimeStep::operator *= (const double x) {
01172     
01173      const unsigned int original_days = _days;      
01174      
01175      const unsigned int original_day_fraction = _day_fraction;      
01176      
01177      std::cerr << "original_days = " << original_days << "  original_day_fraction = " << original_day_fraction << std::endl;
01178      
01179      if (x == 0.0) {
01180      _days = _day_fraction = 0;
01181      _sign = +1;
01182      return * this;
01183      }
01184      
01185      if (x < 0.0) _sign = -_sign;
01186      
01187      _days = (unsigned int)(fabs(x)*_days);
01188      
01189      _day_fraction = (unsigned int)rint((fabs(x)-floor(fabs(x)))*original_day_fraction);
01190      
01191      if (original_days > 0) {
01192      const unsigned int frac = (unsigned int)rint((fabs(x)-floor(fabs(x)))*max_day_fraction());
01193      for (unsigned int k=0;k<original_days;++k) {
01194      _day_fraction += frac;
01195      while (_day_fraction >= max_day_fraction()) {
01196      ++_days;
01197      _day_fraction -= max_day_fraction();
01198      }
01199      }
01200      }
01201      
01202      if (((unsigned int)floor(fabs(x))) > 0) {
01203      // unsigned int count = (unsigned int)floor(fabs(x))-1;
01204      unsigned int count = (unsigned int)floor(fabs(x));
01205      while (count != 0) {
01206      _day_fraction += original_day_fraction;
01207      while (_day_fraction >= max_day_fraction()) {
01208      ++_days;
01209      _day_fraction -= max_day_fraction();
01210      }
01211      --count;
01212      }
01213      }
01214      
01215      std::cerr << "out:   x = " << x << "  _days = " << _days << "  _day_fraction = " << _day_fraction << std::endl;
01216      
01217      return * this;
01218      }
01219   */
01220   
01221   TimeStep TimeStep::operator + () const { return TimeStep(*this); }
01222   
01223   TimeStep TimeStep::operator - () const { return TimeStep(_days,_day_fraction,-_sign); }
01224   
01225   bool TimeStep::operator < (const TimeStep & ts) const {
01226     
01227     if (*this == ts) return false; // they are equal!
01228     
01229     if (_sign == ts._sign) {
01230       if (_sign == -1) {
01231         if (_days > ts._days) return true;
01232         if ( (_days == ts._days) &&
01233              (_day_fraction > ts._day_fraction) ) return true;
01234       } else {
01235         if (_days < ts._days) return true;
01236         if ( (_days == ts._days) &&
01237              (_day_fraction < ts._day_fraction) ) return true;
01238       }
01239     } else {
01240       if (_sign == -1) {
01241         return true;
01242       }
01243     }
01244     return false;
01245   }
01246   
01247   bool TimeStep::operator > (const TimeStep & ts) const {
01248     if (*this == ts) return false; // they are equal!
01249     return (!((*this) < ts));
01250   }
01251   
01252   double TimeStep::GetDouble() const {
01253     return (FromUnits(_sign*(_days+(_day_fraction*1.0)/max_day_fraction()),DAY));
01254   }
01255   
01256   bool TimeStep::operator == (const TimeStep & ts) const {
01257     
01258     // if (ts.IsZero() && IsZero()) return true; // this check is 'probably' not needed since we use the internal_check() call every time we change the object
01259     if (ts._days != _days) return false;
01260     if (ts._day_fraction != _day_fraction) return false;
01261     if (ts._sign != _sign) return false;
01262     return true;
01263   }
01264   
01265   bool TimeStep::operator != (const TimeStep & ts) const {
01266     return (!((*this) == ts));
01267   }
01268   
01269   TimeStep TimeStep::absolute() const { return TimeStep(_days,_day_fraction,+1); }
01270   
01271   bool TimeStep::IsZero() const { return ((_days == 0) && (_day_fraction == 0)); }
01272   
01273   // UniverseTypeAwareTime
01274   
01275   UniverseTypeAwareTime::UniverseTypeAwareTime() { }    
01276   
01277   UniverseTypeAwareTime::UniverseTypeAwareTime(const double t) {
01278     SetTime(t);
01279   }
01280   
01281   UniverseTypeAwareTime::UniverseTypeAwareTime(const Date & d) {
01282     SetTime(d);
01283   }
01284   
01285   UniverseTypeAwareTime::UniverseTypeAwareTime(const UniverseTypeAwareTime & t) {
01286     date = t.date;
01287     time = t.time;
01288   }
01289   
01290   double UniverseTypeAwareTime::GetTime() const {
01291     double _t = 0.0;
01292     switch (universe->GetUniverseType()) {
01293     case Real:
01294       _t = date.GetTime();
01295       break;
01296     case Simulated:
01297       _t = time;
01298       break;
01299     }
01300     return _t;
01301   }
01302   
01303   double UniverseTypeAwareTime::Time() const {
01304     return (GetTime());
01305   }
01306   
01307   Date UniverseTypeAwareTime::GetDate() const {
01308     return (date);
01309   }
01310   
01311   void  UniverseTypeAwareTime::SetTime(const double t) {
01312     date.SetJulian(FromUnits(t,DAY,-1));
01313     time = t;
01314   }
01315   
01316   void UniverseTypeAwareTime::SetTime(const Date & d) {
01317     SetDate(d);
01318   }
01319   
01320   void  UniverseTypeAwareTime::SetDate(const Date & d) {
01321     date = d;
01322     time = d.GetTime();
01323   }
01324   
01325   void UniverseTypeAwareTime::SetTime(const UniverseTypeAwareTime & t) {
01326     date = t.date;
01327     time = t.time;
01328   }
01329   
01330   /* 
01331      UniverseTypeAwareTimeStep UniverseTypeAwareTime::operator + (const UniverseTypeAwareTime & t) const {
01332      switch (universe->GetUniverseType()) {
01333      case Real:
01334      {
01335      UniverseTypeAwareTimeStep _ts(date.GetTimeStep());
01336      _ts += t.GetDate().GetTimeStep();
01337      return _ts;
01338      }
01339      break;
01340      case Simulated:
01341      {
01342      UniverseTypeAwareTimeStep _ts(time);
01343      _ts += t.Time();
01344      return _ts;
01345      }
01346      break;
01347      }
01348      return UniverseTypeAwareTimeStep();
01349      }
01350   */
01351   
01352   UniverseTypeAwareTimeStep UniverseTypeAwareTime::operator - (const UniverseTypeAwareTime & t) const {
01353     switch (universe->GetUniverseType()) {
01354     case Real:
01355       {
01356         UniverseTypeAwareTimeStep _ts(date.GetTimeStep());
01357         _ts -= t.GetDate().GetTimeStep();
01358         return _ts;
01359       }
01360       break;
01361     case Simulated:
01362       {
01363         UniverseTypeAwareTimeStep _ts(time);
01364         _ts -= t.Time();
01365         return _ts;
01366       }
01367       break;
01368     }
01369     return UniverseTypeAwareTimeStep();
01370   }
01371   
01372   /* 
01373      UniverseTypeAwareTimeStep UniverseTypeAwareTime::operator + (const UniverseTypeAwareTimeStep & ts) const {
01374      switch (universe->GetUniverseType()) {
01375      case Real:
01376      {
01377      UniverseTypeAwareTimeStep _ts(date.GetTimeStep());
01378      _ts += ts;
01379      return _ts;
01380      }
01381      break;
01382      case Simulated:
01383      {
01384      UniverseTypeAwareTimeStep _ts(time);
01385      _ts += ts;
01386      return _ts;
01387      }
01388      break;
01389      }
01390      return UniverseTypeAwareTimeStep();
01391      }
01392   */
01393   
01394   /* 
01395      UniverseTypeAwareTimeStep UniverseTypeAwareTime::operator - (const UniverseTypeAwareTimeStep & ts) const {
01396      switch (universe->GetUniverseType()) {
01397      case Real:
01398      {
01399      UniverseTypeAwareTimeStep _ts(date.GetTimeStep());
01400      _ts -= ts;
01401      return _ts;
01402      }
01403      break;
01404      case Simulated:
01405      {
01406      UniverseTypeAwareTimeStep _ts(time);
01407      _ts -= ts;
01408      return _ts;
01409      }
01410      break;
01411      }
01412      return UniverseTypeAwareTimeStep();
01413      }
01414   */
01415   
01416   UniverseTypeAwareTime UniverseTypeAwareTime::operator + (const UniverseTypeAwareTimeStep & ts) const {
01417     switch (universe->GetUniverseType()) {
01418     case Real:
01419       {
01420         UniverseTypeAwareTime _t(date);
01421         _t += ts;
01422         return _t;
01423       }
01424       break;
01425     case Simulated:
01426       {
01427         UniverseTypeAwareTime _t(time);
01428         _t += ts;
01429         return _t;
01430       }
01431       break;
01432     }
01433     return UniverseTypeAwareTime();
01434   }
01435   
01436   UniverseTypeAwareTime UniverseTypeAwareTime::operator - (const UniverseTypeAwareTimeStep & ts) const {
01437     switch (universe->GetUniverseType()) {
01438     case Real:
01439       {
01440         UniverseTypeAwareTime _t(date);
01441         _t -= ts;
01442         return _t;
01443       }
01444       break;
01445     case Simulated:
01446       {
01447         UniverseTypeAwareTime _t(time);
01448         _t -= ts;
01449         return _t;
01450       }
01451       break;
01452     }
01453     return UniverseTypeAwareTime();
01454   }
01455   
01456   UniverseTypeAwareTime & UniverseTypeAwareTime::operator += (const UniverseTypeAwareTimeStep & ts_in) {
01457     switch (universe->GetUniverseType()) {
01458     case Real:
01459       date += ts_in;
01460       break;
01461     case Simulated:
01462       time += ts_in.GetDouble();
01463       break;
01464     }
01465     return * this;
01466   }
01467   
01468   UniverseTypeAwareTime & UniverseTypeAwareTime::operator -= (const UniverseTypeAwareTimeStep & ts_in) {
01469     switch (universe->GetUniverseType()) {
01470     case Real:
01471       date -= ts_in;
01472       break;
01473     case Simulated:
01474       time -= ts_in.GetDouble();
01475       break;
01476     }
01477     return * this;
01478   }
01479   
01480   bool UniverseTypeAwareTime::operator == (const UniverseTypeAwareTime & t) const {
01481     bool _z = false;
01482     switch (universe->GetUniverseType()) {
01483     case Real:
01484       // if (date.GetTimeStep() == t.GetDate().GetTimeStep()) _z = true;
01485       if (date == t.GetDate()) _z = true;
01486       break;
01487     case Simulated:
01488       if (time == t.time) _z = true;
01489       break;
01490     }
01491     return _z;
01492   }
01493   
01494   bool UniverseTypeAwareTime::operator != (const UniverseTypeAwareTime & t) const {
01495     return (!((*this) == t));
01496   }
01497   
01498   bool UniverseTypeAwareTime::operator < (const UniverseTypeAwareTime & t) const {
01499     if (*this == t) return false;
01500     bool _z = false;
01501     switch (universe->GetUniverseType()) {
01502     case Real:
01503       if (date.GetTimeStep() < t.GetDate().GetTimeStep()) _z = true;
01504       break;
01505     case Simulated:
01506       if (time < t.time) _z = true;
01507       break;
01508     }
01509     return _z;
01510   }
01511   
01512   bool UniverseTypeAwareTime::operator > (const UniverseTypeAwareTime & t) const {
01513     if (*this == t) return false;
01514     bool _z = false;
01515     switch (universe->GetUniverseType()) {
01516     case Real:
01517       if (date.GetTimeStep() > t.GetDate().GetTimeStep()) _z = true;
01518       break;
01519     case Simulated:
01520       if (time > t.time) _z = true;
01521       break;
01522     }
01523     return _z;
01524   }
01525   
01526   bool UniverseTypeAwareTime::operator <= (const UniverseTypeAwareTime & t) const {
01527     bool _z = false;
01528     switch (universe->GetUniverseType()) {
01529     case Real:
01530       if (date.GetTimeStep() <= t.GetDate().GetTimeStep()) _z = true;
01531       break;
01532     case Simulated:
01533       if (time <= t.time) _z = true;
01534       break;
01535     }
01536     return _z;
01537   }
01538   
01539   bool UniverseTypeAwareTime::operator >= (const UniverseTypeAwareTime & t) const {
01540     bool _z = false;
01541     switch (universe->GetUniverseType()) {
01542     case Real:
01543       if (date.GetTimeStep() >= t.GetDate().GetTimeStep()) _z = true;
01544       break;
01545     case Simulated:
01546       if (time >= t.time) _z = true;
01547       break;
01548     }
01549     return _z;
01550   }
01551   
01552   // out...
01553   
01554   UniverseTypeAwareTimeStep operator * (const double x, const UniverseTypeAwareTimeStep & ts) {
01555     UniverseTypeAwareTimeStep _ts(ts);
01556     _ts *= x;
01557     return _ts;
01558   }
01559   
01560   UniverseTypeAwareTimeStep operator * (const UniverseTypeAwareTimeStep & ts, const double x) {
01561     UniverseTypeAwareTimeStep _ts(ts);
01562     _ts *= x;
01563     return _ts; 
01564   }
01565   
01566   UniverseTypeAwareTimeStep operator * (const int i, const UniverseTypeAwareTimeStep & ts) {
01567     UniverseTypeAwareTimeStep _ts(ts);
01568     _ts *= i;
01569     return _ts; 
01570   }
01571   
01572   UniverseTypeAwareTimeStep operator * (const UniverseTypeAwareTimeStep & ts, const int i) {
01573     UniverseTypeAwareTimeStep _ts(ts);
01574     _ts *= i;
01575     return _ts; 
01576   }
01577   
01578   // UniverseTypeAwareTimeStep 
01579 
01580   UniverseTypeAwareTimeStep::UniverseTypeAwareTimeStep() : ts(0,0,1), dts(0.0) {
01581     
01582   }
01583   
01584   UniverseTypeAwareTimeStep::UniverseTypeAwareTimeStep(const double t) : ts(t), dts(t) {
01585     
01586   }
01587   
01588   UniverseTypeAwareTimeStep::UniverseTypeAwareTimeStep(const TimeStep & ts_in) : ts(ts_in), dts(ts_in.GetDouble()) {
01589     
01590   }
01591   
01592   UniverseTypeAwareTimeStep::UniverseTypeAwareTimeStep(const int days, const unsigned int day_fraction, const int sign) : ts(days,day_fraction,sign), dts(ts.GetDouble()) {
01593     
01594   }
01595   
01596   UniverseTypeAwareTimeStep::UniverseTypeAwareTimeStep(const UniverseTypeAwareTimeStep & ts_in) : ts(ts_in.ts), dts(ts_in.dts) {
01597     
01598   }
01599   
01600   UniverseTypeAwareTimeStep &  UniverseTypeAwareTimeStep::operator += (const UniverseTypeAwareTimeStep & ts_in) {
01601     ts  += ts_in.ts;
01602     dts += ts_in.dts;
01603     return * this;
01604   }
01605   
01606   UniverseTypeAwareTimeStep &  UniverseTypeAwareTimeStep::operator -= (const UniverseTypeAwareTimeStep & ts_in) {
01607     ts  -= ts_in.ts;
01608     dts -= ts_in.dts;
01609     return * this;
01610   }
01611   
01612   UniverseTypeAwareTimeStep & UniverseTypeAwareTimeStep::operator *= (const int p) {
01613     ts  *= p;
01614     dts *= p;
01615     return * this;
01616   }
01617   
01618   UniverseTypeAwareTimeStep & UniverseTypeAwareTimeStep::operator *= (const double x) {
01619     ts  *= x;
01620     dts *= x;
01621     return * this;
01622   }     
01623   
01624   double UniverseTypeAwareTimeStep::GetDouble() const {
01625     double _x = 0.0;
01626     switch (universe->GetUniverseType()) {
01627     case Real:
01628       _x = ts.GetDouble();
01629       break;
01630     case Simulated:
01631       _x = dts;
01632       break;
01633     }
01634     return _x;
01635   }
01636   
01637   bool UniverseTypeAwareTimeStep::operator < (const double x) const {
01638     return (GetDouble() < x);
01639   }
01640   
01641   bool UniverseTypeAwareTimeStep::operator > (const double x) const {
01642     return (GetDouble() > x);
01643   }
01644   
01645   bool UniverseTypeAwareTimeStep::operator < (const UniverseTypeAwareTime & t) const {
01646     bool _z = false;
01647     switch (universe->GetUniverseType()) {
01648     case Real:
01649       if (ts < t.GetDate().GetTimeStep()) _z = true;
01650       break;
01651     case Simulated:
01652       if (dts < t.GetTime()) _z = true;
01653       break;
01654     }
01655     return _z;
01656   }
01657   
01658   bool UniverseTypeAwareTimeStep::operator > (const UniverseTypeAwareTime & t) const {
01659     bool _z = false;
01660     switch (universe->GetUniverseType()) {
01661     case Real:
01662       if (ts > t.GetDate().GetTimeStep()) _z = true;
01663       break;
01664     case Simulated:
01665       if (dts > t.GetTime()) _z = true;
01666       break;
01667     }
01668     return _z;
01669   }
01670   
01671   UniverseTypeAwareTimeStep UniverseTypeAwareTimeStep::operator + () const {
01672     return * this;
01673   }
01674   
01675   UniverseTypeAwareTimeStep UniverseTypeAwareTimeStep::operator - () const {
01676     UniverseTypeAwareTimeStep _ts;
01677     _ts.ts  = - ts;
01678     _ts.dts = - dts;
01679     return _ts;
01680   }
01681   
01682   UniverseTypeAwareTimeStep UniverseTypeAwareTimeStep::operator + (const UniverseTypeAwareTimeStep & ts_in) const {
01683     UniverseTypeAwareTimeStep _ts(*this);
01684     _ts.ts  += ts_in.ts;
01685     _ts.dts += ts_in.dts;
01686     return _ts;
01687   }
01688   
01689   UniverseTypeAwareTimeStep UniverseTypeAwareTimeStep::operator - (const UniverseTypeAwareTimeStep & ts_in) const {
01690     UniverseTypeAwareTimeStep _ts(*this);
01691     _ts.ts  -= ts_in.ts;
01692     _ts.dts -= ts_in.dts;
01693     return _ts;
01694   }
01695   
01696   UniverseTypeAwareTimeStep UniverseTypeAwareTimeStep::operator + (const UniverseTypeAwareTime & t) const {
01697     UniverseTypeAwareTimeStep _ts(*this);
01698     switch (universe->GetUniverseType()) {
01699     case Real:
01700       _ts += t.GetDate().GetTimeStep();
01701       break;
01702     case Simulated:
01703       _ts.dts += t.GetTime();
01704       break;
01705     }
01706     return _ts;
01707   }
01708   
01709   UniverseTypeAwareTimeStep UniverseTypeAwareTimeStep::operator - (const UniverseTypeAwareTime & t) const {
01710     UniverseTypeAwareTimeStep _ts(*this);
01711     switch (universe->GetUniverseType()) {
01712     case Real:
01713       _ts.ts -= t.GetDate().GetTimeStep();
01714       break;
01715     case Simulated:
01716       _ts.dts -= t.GetTime();
01717       break;
01718     }
01719     return _ts;
01720   }
01721   
01722   UniverseTypeAwareTimeStep UniverseTypeAwareTimeStep::absolute() const {
01723     UniverseTypeAwareTimeStep _ts(*this);
01724     switch (universe->GetUniverseType()) {
01725     case Real:
01726       _ts.ts = ts.absolute();
01727       break;
01728     case Simulated:
01729       _ts.dts = fabs(dts);
01730       break;
01731     }
01732     return _ts;
01733   }
01734   
01735   bool UniverseTypeAwareTimeStep::IsZero() const {
01736     bool _z = false;
01737     switch (universe->GetUniverseType()) {
01738     case Real:
01739       _z = ts.IsZero();
01740       break;
01741     case Simulated:
01742       _z = (dts == 0.0);
01743       break;
01744     }
01745     return _z;
01746   }
01747   
01748   bool UniverseTypeAwareTimeStep::operator < (const UniverseTypeAwareTimeStep & ts_in) const {
01749     if (*this == ts_in) return false;
01750     bool _z = false;
01751     switch (universe->GetUniverseType()) {
01752     case Real:
01753       if (ts < ts_in.ts) _z = true;
01754       break;
01755     case Simulated:
01756       if (dts < ts_in.dts) _z = true;
01757       break;
01758     }
01759     return _z;
01760   }
01761   
01762   bool UniverseTypeAwareTimeStep::operator > (const UniverseTypeAwareTimeStep & ts_in) const {
01763     if (*this == ts_in) return false;
01764     bool _z = false;
01765     switch (universe->GetUniverseType()) {
01766     case Real:
01767       if (ts > ts_in.ts) _z = true;
01768       break;
01769     case Simulated:
01770       if (dts > ts_in.dts) _z = true;
01771       break;
01772     }
01773     return _z; 
01774   }
01775   
01776   bool UniverseTypeAwareTimeStep::operator == (const UniverseTypeAwareTimeStep & ts_in) const {
01777     bool _z = false;
01778     switch (universe->GetUniverseType()) {
01779     case Real:
01780       if (ts == ts_in.ts) _z = true;
01781       break;
01782     case Simulated:
01783       if (dts == ts_in.dts) _z = true;
01784       break;
01785     }
01786     return _z;
01787   }
01788   
01789   bool UniverseTypeAwareTimeStep::operator != (const UniverseTypeAwareTimeStep & ts_in) const {
01790     return (!((*this) == ts_in));
01791   }
01792   
01793   // Angle
01794   
01795   void Angle::SetRad(double a) {
01796     radians = a;
01797   }
01798   
01799   void Angle::GetRad(double & a) const {
01800     a = radians;
01801   }
01802   
01803   double Angle::GetRad() const {
01804     return radians;
01805   }
01806   
01807   void Angle::SetDPS(double  d, double  p, double  s) {
01808     if (d >= 0) 
01809       radians = (pi/180)*(d+p/60.0+s/3600.0);
01810     else 
01811       radians = (pi/180)*(d-p/60.0-s/3600.0);
01812   }
01813   
01814   void Angle::GetDPS(double &d, double &p, double &s) const {
01815     double frac, fdeg = (180/pi)*radians;
01816     if (fdeg < 0.0) {
01817       d    = -floor(-fdeg);
01818       frac = d - fdeg;
01819     } else {
01820       d    = floor(fdeg);
01821       frac = fdeg - d;
01822     }
01823     //
01824     p = floor(frac*60.0);  
01825     s = frac*3600.0 - p*60.0;
01826   } 
01827   
01828   void Angle::SetHMS(double  h, double  m, double  s) {
01829     if (h >= 0) 
01830       radians = 15*(pi/180)*(h+m/60.0+s/3600.0);
01831     else
01832       radians = 15*(pi/180)*(h-m/60.0-s/3600.0);
01833   }
01834   
01835   void Angle::GetHMS(double &h, double &m, double &s) const {
01836     double frac, fh = (180/pi)*radians/15.0;
01837     if (fh < 0.0) {
01838       h    = -floor(-fh);
01839       frac = h - fh;
01840     } else {
01841       h    = floor(fh);
01842       frac = fh - h;
01843     }
01844     //
01845     m = floor(frac*60.0);  
01846     s = frac*3600.0 - m*60.0; 
01847   }    
01848   
01849   // reference systems
01850   
01851   Angle obleq(const Date &date) {
01852     // T in centuries from J2000 
01853     Date d = date;
01854     // d.ConvertToTimeScale(UT); // UT or UTC ??
01855     // const double T = (d.GetJulian() - 2451545.0)/36525.0;
01856     // DOUBLE-CHECK this "UT"!!!
01857     const double T = (d.GetJulian(UT) - 2451545.0)/36525.0;
01858     Angle a;
01859     // updated Feb 2004
01860     a.SetDPS(23,26,21.448+((0.001813*T-0.00059)*T-46.8150)*T);
01861     return a;
01862   }
01863   
01864   // short term correction
01865   /* 
01866      Angle delta_obleq(const Date &date) {
01867      Date d = date;
01868      d.ConvertToTimeScale(UT); // UT or UTC ??
01869      // days from 2451545.0
01870      const double days = (d.GetJulian() - 2451545.0);
01871      Angle a;
01872      // updated Feb 2004
01873      const double delta_grad = ( 0.0026*cos((pi/180.0)*(125.0-0.05295*days)) +
01874      0.0002*cos((pi/180.0)*(200.9+1.97129*days)) );
01875      a.SetDPS(delta_grad,0,0);
01876      return a;
01877      }
01878   */
01879   
01880   // short term longitude correction
01881   /* 
01882      Angle delta_longitude(const Date &date) {
01883      Date d = date;
01884      d.ConvertToTimeScale(UT); // UT or UTC ??
01885      // days from 2451545.0
01886      const double days = (d.GetJulian() - 2451545.0);
01887      Angle a;
01888      // updated Feb 2004
01889      const double delta_grad = ( -0.0048*sin((pi/180.0)*(125.0-0.05295*days))
01890      -0.0004*sin((pi/180.0)*(200.9+1.97129*days)) );
01891      a.SetDPS(delta_grad,0,0);
01892      return a;
01893      }
01894   */
01895   
01896   Angle gmst(const Date &date) {
01897     // tested using xephem, very good agreement!
01898     // T in centuries from JD 2451545.0 UT1=UT
01899     Date d = date;
01900     // d.ConvertToTimeScale(UT); // UT or UTC ??
01901     // const double T = (d.GetJulian() - 2451545.0)/36525.0;
01902     const double T = (d.GetJulian(UT) - 2451545.0)/36525.0;
01903     Angle a;
01904     a.SetHMS(d.GetDayFraction()*24+6,41,50.54841+((-0.0000062*T+0.093104)*T+8640184.812866)*T);
01905     return a;
01906   }
01907   
01908   // J2000 versions
01909   Angle obleq_J2000() {
01910     Date d; d.SetJ2000();
01911     return obleq(d);
01912   }
01913   
01914   /* 
01915      Angle delta_obleq_J2000() {
01916      Date d; d.SetJ2000();
01917      return delta_obleq(d);
01918      }
01919      
01920      Angle delta_longitude_J2000() {
01921      Date d; d.SetJ2000();
01922      return delta_longitude(d);
01923      }
01924   */
01925   
01926   //
01927   
01928   /* 
01929      void EquatorialToDate_To_EquatorialJ2000(Vector &v, const Date &date) {
01930      // to be checked!
01931      }
01932   */
01933   
01934   //
01935   
01936   /* 
01937      void EclipticToEquatorial(Vector &v, const Date &date) {
01938      v.rotate(0,0,delta_longitude(date).GetRad());
01939      v.rotate(0,obleq(date).GetRad()+delta_obleq(date).GetRad(),0);  
01940      }
01941   */
01942   //
01943   /* 
01944      void EclipticToEquatorial(Vector &v, const Date &date) {
01945      v.rotate(0,0,delta_longitude(date).GetRad());
01946      v.rotate(0,obleq(date).GetRad()+delta_obleq(date).GetRad(),0);  
01947      }
01948      
01949      void EquatorialToEcliptic(Vector &v, const Date &date) {
01950      v.rotate(0,-obleq(date).GetRad()-delta_obleq(date).GetRad(),0);  
01951      v.rotate(0,0,-delta_longitude(date).GetRad());
01952      }
01953      
01954      void EclipticToEquatorial_J2000(Vector &v) {
01955      v.rotate(0,0,delta_longitude_J2000().GetRad());
01956      v.rotate(0,obleq_J2000().GetRad()+delta_obleq_J2000().GetRad(),0);  
01957      }
01958      
01959      void EquatorialToEcliptic_J2000(Vector &v) {
01960      v.rotate(0,-obleq_J2000().GetRad()-delta_obleq_J2000().GetRad(),0);  
01961      v.rotate(0,0,-delta_longitude_J2000().GetRad());
01962      }
01963   */
01964   
01965   void EclipticToEquatorial(Vector &v, const Date &date) {
01966     v.rotate(0,obleq(date).GetRad(),0);  
01967   }
01968   
01969   void EquatorialToEcliptic(Vector &v, const Date &date) {
01970     v.rotate(0,-obleq(date).GetRad(),0);  
01971   }
01972   
01973   void EclipticToEquatorial_J2000(Vector &v) {
01974     v.rotate(0,obleq_J2000().GetRad(),0);  
01975   }
01976   
01977   void EquatorialToEcliptic_J2000(Vector &v) {
01978     v.rotate(0,-obleq_J2000().GetRad(),0);  
01979   }
01980   
01981   // data from Table II, same paper of alpha_delta_meridian(...);
01982   static void alpha_delta_meridian_moon(const Date & date,
01983                                         Angle & alpha_zero, Angle & delta_zero, Angle & W) {
01984     Date tmp_date; 
01985     tmp_date.SetJ2000();
01986     const Date date_J2000(tmp_date);
01987     
01988     const double d = (date.GetJulian()-date_J2000.GetJulian());
01989     const double T = d / 36525.0;
01990     
01991     Angle E1, E2, E3, E4, E5, E6, E7, E8, E9, E10, E11, E12, E13;
01992     
01993     E1.SetDPS( 125.045 -  0.0529921*d,0,0);
01994     E2.SetDPS( 250.089 -  0.1059842*d,0,0);
01995     E3.SetDPS( 260.008 + 13.0120009*d,0,0);
01996     E4.SetDPS( 176.625 + 13.3407154*d,0,0);
01997     E5.SetDPS( 357.529 +  0.9856003*d,0,0);
01998     E6.SetDPS( 311.589 + 26.4057084*d,0,0);
01999     E7.SetDPS( 134.963 + 13.0649930*d,0,0);
02000     E8.SetDPS( 276.617 +  0.3287146*d,0,0);
02001     E9.SetDPS(  34.226 +  1.7484877*d,0,0);
02002     E10.SetDPS( 15.134 -  0.1589763*d,0,0);
02003     E11.SetDPS(119.743 +  0.0036096*d,0,0);
02004     E12.SetDPS(239.961 +  0.1643573*d,0,0);
02005     E13.SetDPS( 25.053 + 12.9590088*d,0,0);
02006     
02007     alpha_zero.SetDPS(269.9949 
02008                       + 0.0031*T
02009                       - 3.8787*sin(E1)
02010                       - 0.1204*sin(E2)
02011                       + 0.0700*sin(E3)
02012                       - 0.0172*sin(E4)
02013                       + 0.0072*sin(E6)
02014                       - 0.0052*sin(E10)
02015                       + 0.0043*sin(E13),0,0);
02016     
02017     delta_zero.SetDPS(66.5392
02018                       + 0.0130*T
02019                       + 1.5419*cos(E1)
02020                       + 0.0239*cos(E2)
02021                       - 0.0278*cos(E3)
02022                       + 0.0068*cos(E4)
02023                       - 0.0029*cos(E6)
02024                       + 0.0009*cos(E7)
02025                       + 0.0008*cos(E10)
02026                       - 0.0009*cos(E13),0,0);
02027     
02028     W.SetDPS(38.3213
02029              + 13.17635815*d
02030              -  1.4e-12*d*d
02031              +  3.5610*sin(E1)
02032              +  0.1208*sin(E2)
02033              -  0.0642*sin(E3)
02034              +  0.0158*sin(E4)
02035              +  0.0252*sin(E5)
02036              -  0.0066*sin(E6)
02037              -  0.0047*sin(E7)
02038              -  0.0046*sin(E8)
02039              +  0.0028*sin(E9)
02040              +  0.0052*sin(E10)
02041              +  0.0040*sin(E11)
02042              +  0.0019*sin(E12)
02043              -  0.0044*sin(E13),0,0);
02044   }
02045   
02046   void alpha_delta_meridian(const JPL_planets p, const Date & date,
02047                             Angle & alpha_zero, Angle & delta_zero, Angle & W) {
02048     
02049     Date tmp_date; 
02050     tmp_date.SetJ2000();
02051     const Date date_J2000(tmp_date);
02052     
02053     // CHECK TIMESCALES!!
02054     const double d = (date.GetJulian()-date_J2000.GetJulian());
02055     const double T = d / 36525.0;
02056     
02057     switch (p) {
02058     case SUN:   
02059       alpha_zero.SetDPS(286.13,0,0);
02060       delta_zero.SetDPS(63.87,0,0);
02061       W.SetDPS(84.10+14.1844000*d,0,0);
02062       break;
02063     case MERCURY: 
02064       alpha_zero.SetDPS(281.01-0.033*T,0,0);
02065       delta_zero.SetDPS(61.45-0.005*T,0,0);
02066       W.SetDPS(329.548+6.1385025*d,0,0);
02067       break;
02068     case VENUS:  
02069       alpha_zero.SetDPS(272.76,0,0);
02070       delta_zero.SetDPS(67.16,0,0);
02071       W.SetDPS(160.20-1.4813688*d,0,0);
02072       break;
02073     case EARTH:
02074       alpha_zero.SetDPS(0.0-0.641*T,0,0);
02075       delta_zero.SetDPS(90.0-0.557*T,0,0);
02076       W.SetDPS(190.147+360.9856235*d,0,0);
02077       break;
02078     case MOON:
02079       alpha_delta_meridian_moon(date,alpha_zero,delta_zero,W);
02080       break;
02081     case MARS:   
02082       alpha_zero.SetDPS(317.68143-0.1061*T,0,0);
02083       delta_zero.SetDPS(52.88650-0.0609*T,0,0);
02084       W.SetDPS(176.630+350.89198226*d,0,0);
02085       break;
02086     case JUPITER: 
02087       alpha_zero.SetDPS(268.05-0.009*T,0,0);
02088       delta_zero.SetDPS(64.49+0.003*T,0,0);
02089       W.SetDPS(284.95+870.5366420*d,0,0);
02090       break;
02091     case SATURN:  
02092       alpha_zero.SetDPS(40.589-0.036*T,0,0);
02093       delta_zero.SetDPS(83.537-0.004*T,0,0);
02094       W.SetDPS(38.90+810.7939024*d,0,0);
02095       break;
02096     case URANUS: 
02097       alpha_zero.SetDPS(257.311,0,0);
02098       delta_zero.SetDPS(-15.175,0,0);
02099       W.SetDPS(203.81-501.1600928*d,0,0);
02100       break;
02101     case NEPTUNE:
02102       {
02103         Angle N;
02104         N.SetDPS(357.85+52.316*T,0,0);
02105         //
02106         alpha_zero.SetDPS(299.36+0.70*sin(N),0,0);
02107         delta_zero.SetDPS(43.46-0.51*cos(N),0,0);
02108         W.SetDPS(253.18+536.3128492*d-0.48*sin(N),0,0);
02109       }
02110       break;
02111     case PLUTO: 
02112       alpha_zero.SetDPS(313.02,0,0);
02113       delta_zero.SetDPS(9.09,0,0);
02114       W.SetDPS(236.77-56.3623195*d,0,0);
02115       break;
02116       //
02117     default: 
02118       alpha_zero.SetRad(0.0);
02119       delta_zero.SetRad(0.0);
02120       W.SetRad(0.0);
02121       break;
02122     }
02123   }
02124   
02125 } // namespace orsa
02126 

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