00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "orsa_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
00098
00099
00100
00101
00102 if (config->paths[JPL_EPHEM_FILE]->GetValue() != "") {
00103
00104
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
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";
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
00216
00217
00218 TimeScale default_Date_timescale = TT;
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
00239
00240
00241
00242
00243
00244
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},
00272 {0,0,0,0}
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
00294
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}
00495 };
00496
00497
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
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
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
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
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
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
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
00722
00723
00724
00725
00726
00727
00728
00729
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));
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));
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
00870
00871
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,
00879 UTC);
00880 }
00881
00882 void Date::SetToday() {
00883 SetNow();
00884 df = 0;
00885 }
00886
00887 void Date::SetJ2000() {
00888 SetJulian(2451545.0,TT);
00889
00890
00891
00892
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
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
00989
00990
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
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
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
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
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134 TimeStep & TimeStep::operator *= (const int p) {
01135 const double t = GetDouble()*p;
01136
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
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
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
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;
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;
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
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
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
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
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
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
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
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
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
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
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
01850
01851 Angle obleq(const Date &date) {
01852
01853 Date d = date;
01854
01855
01856
01857 const double T = (d.GetJulian(UT) - 2451545.0)/36525.0;
01858 Angle a;
01859
01860 a.SetDPS(23,26,21.448+((0.001813*T-0.00059)*T-46.8150)*T);
01861 return a;
01862 }
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896 Angle gmst(const Date &date) {
01897
01898
01899 Date d = date;
01900
01901
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
01909 Angle obleq_J2000() {
01910 Date d; d.SetJ2000();
01911 return obleq(d);
01912 }
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
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
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
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 }
02126