1 /************************************************************************ 2 * Copyright (C) 1996-2008, International Business Machines Corporation * 3 * and others. All Rights Reserved. * 4 ************************************************************************ 5 * 2003-nov-07 srl Port from Java 6 */ 7 8 #include "astro.h" 9 10 #if !UCONFIG_NO_FORMATTING 11 12 #include "unicode/calendar.h" 13 #include <math.h> 14 #include <float.h> 15 #include "unicode/putil.h" 16 #include "uhash.h" 17 #include "umutex.h" 18 #include "ucln_in.h" 19 #include "putilimp.h" 20 #include <stdio.h> // for toString() 21 22 #if defined (PI) 23 #undef PI 24 #endif 25 26 #ifdef U_DEBUG_ASTRO 27 # include "uresimp.h" // for debugging 28 29 static void debug_astro_loc(const char *f, int32_t l) 30 { 31 fprintf(stderr, "%s:%d: ", f, l); 32 } 33 34 static void debug_astro_msg(const char *pat, ...) 35 { 36 va_list ap; 37 va_start(ap, pat); 38 vfprintf(stderr, pat, ap); 39 fflush(stderr); 40 } 41 #include "unicode/datefmt.h" 42 #include "unicode/ustring.h" 43 static const char * debug_astro_date(UDate d) { 44 static char gStrBuf[1024]; 45 static DateFormat *df = NULL; 46 if(df == NULL) { 47 df = DateFormat::createDateTimeInstance(DateFormat::MEDIUM, DateFormat::MEDIUM, Locale::getUS()); 48 df->adoptTimeZone(TimeZone::getGMT()->clone()); 49 } 50 UnicodeString str; 51 df->format(d,str); 52 u_austrncpy(gStrBuf,str.getTerminatedBuffer(),sizeof(gStrBuf)-1); 53 return gStrBuf; 54 } 55 56 // must use double parens, i.e.: U_DEBUG_ASTRO_MSG(("four is: %d",4)); 57 #define U_DEBUG_ASTRO_MSG(x) {debug_astro_loc(__FILE__,__LINE__);debug_astro_msg x;} 58 #else 59 #define U_DEBUG_ASTRO_MSG(x) 60 #endif 61 62 static inline UBool isINVALID(double d) { 63 return(uprv_isNaN(d)); 64 } 65 66 static UMTX ccLock = NULL; 67 68 U_CDECL_BEGIN 69 static UBool calendar_astro_cleanup(void) { 70 umtx_destroy(&ccLock); 71 return TRUE; 72 } 73 U_CDECL_END 74 75 U_NAMESPACE_BEGIN 76 77 /** 78 * The number of standard hours in one sidereal day. 79 * Approximately 24.93. 80 * @internal 81 * @deprecated ICU 2.4. This class may be removed or modified. 82 */ 83 #define SIDEREAL_DAY (23.93446960027) 84 85 /** 86 * The number of sidereal hours in one mean solar day. 87 * Approximately 24.07. 88 * @internal 89 * @deprecated ICU 2.4. This class may be removed or modified. 90 */ 91 #define SOLAR_DAY (24.065709816) 92 93 /** 94 * The average number of solar days from one new moon to the next. This is the time 95 * it takes for the moon to return the same ecliptic longitude as the sun. 96 * It is longer than the sidereal month because the sun's longitude increases 97 * during the year due to the revolution of the earth around the sun. 98 * Approximately 29.53. 99 * 100 * @see #SIDEREAL_MONTH 101 * @internal 102 * @deprecated ICU 2.4. This class may be removed or modified. 103 */ 104 const double CalendarAstronomer::SYNODIC_MONTH = 29.530588853; 105 106 /** 107 * The average number of days it takes 108 * for the moon to return to the same ecliptic longitude relative to the 109 * stellar background. This is referred to as the sidereal month. 110 * It is shorter than the synodic month due to 111 * the revolution of the earth around the sun. 112 * Approximately 27.32. 113 * 114 * @see #SYNODIC_MONTH 115 * @internal 116 * @deprecated ICU 2.4. This class may be removed or modified. 117 */ 118 #define SIDEREAL_MONTH 27.32166 119 120 /** 121 * The average number number of days between successive vernal equinoxes. 122 * Due to the precession of the earth's 123 * axis, this is not precisely the same as the sidereal year. 124 * Approximately 365.24 125 * 126 * @see #SIDEREAL_YEAR 127 * @internal 128 * @deprecated ICU 2.4. This class may be removed or modified. 129 */ 130 #define TROPICAL_YEAR 365.242191 131 132 /** 133 * The average number of days it takes 134 * for the sun to return to the same position against the fixed stellar 135 * background. This is the duration of one orbit of the earth about the sun 136 * as it would appear to an outside observer. 137 * Due to the precession of the earth's 138 * axis, this is not precisely the same as the tropical year. 139 * Approximately 365.25. 140 * 141 * @see #TROPICAL_YEAR 142 * @internal 143 * @deprecated ICU 2.4. This class may be removed or modified. 144 */ 145 #define SIDEREAL_YEAR 365.25636 146 147 //------------------------------------------------------------------------- 148 // Time-related constants 149 //------------------------------------------------------------------------- 150 151 /** 152 * The number of milliseconds in one second. 153 * @internal 154 * @deprecated ICU 2.4. This class may be removed or modified. 155 */ 156 #define SECOND_MS U_MILLIS_PER_SECOND 157 158 /** 159 * The number of milliseconds in one minute. 160 * @internal 161 * @deprecated ICU 2.4. This class may be removed or modified. 162 */ 163 #define MINUTE_MS U_MILLIS_PER_MINUTE 164 165 /** 166 * The number of milliseconds in one hour. 167 * @internal 168 * @deprecated ICU 2.4. This class may be removed or modified. 169 */ 170 #define HOUR_MS U_MILLIS_PER_HOUR 171 172 /** 173 * The number of milliseconds in one day. 174 * @internal 175 * @deprecated ICU 2.4. This class may be removed or modified. 176 */ 177 #define DAY_MS U_MILLIS_PER_DAY 178 179 /** 180 * The start of the julian day numbering scheme used by astronomers, which 181 * is 1/1/4713 BC (Julian), 12:00 GMT. This is given as the number of milliseconds 182 * since 1/1/1970 AD (Gregorian), a negative number. 183 * Note that julian day numbers and 184 * the Julian calendar are <em>not</em> the same thing. Also note that 185 * julian days start at <em>noon</em>, not midnight. 186 * @internal 187 * @deprecated ICU 2.4. This class may be removed or modified. 188 */ 189 #define JULIAN_EPOCH_MS -210866760000000.0 190 191 192 /** 193 * Milliseconds value for 0.0 January 2000 AD. 194 */ 195 #define EPOCH_2000_MS 946598400000.0 196 197 //------------------------------------------------------------------------- 198 // Assorted private data used for conversions 199 //------------------------------------------------------------------------- 200 201 // My own copies of these so compilers are more likely to optimize them away 202 const double CalendarAstronomer::PI = 3.14159265358979323846; 203 204 #define CalendarAstronomer_PI2 (CalendarAstronomer::PI*2.0) 205 #define RAD_HOUR ( 12 / CalendarAstronomer::PI ) // radians -> hours 206 #define DEG_RAD ( CalendarAstronomer::PI / 180 ) // degrees -> radians 207 #define RAD_DEG ( 180 / CalendarAstronomer::PI ) // radians -> degrees 208 209 /*** 210 * Given 'value', add or subtract 'range' until 0 <= 'value' < range. 211 * The modulus operator. 212 */ 213 inline static double normalize(double value, double range) { 214 return value - range * ClockMath::floorDivide(value, range); 215 } 216 217 /** 218 * Normalize an angle so that it's in the range 0 - 2pi. 219 * For positive angles this is just (angle % 2pi), but the Java 220 * mod operator doesn't work that way for negative numbers.... 221 */ 222 inline static double norm2PI(double angle) { 223 return normalize(angle, CalendarAstronomer::PI * 2.0); 224 } 225 226 /** 227 * Normalize an angle into the range -PI - PI 228 */ 229 inline static double normPI(double angle) { 230 return normalize(angle + CalendarAstronomer::PI, CalendarAstronomer::PI * 2.0) - CalendarAstronomer::PI; 231 } 232 233 //------------------------------------------------------------------------- 234 // Constructors 235 //------------------------------------------------------------------------- 236 237 /** 238 * Construct a new <code>CalendarAstronomer</code> object that is initialized to 239 * the current date and time. 240 * @internal 241 * @deprecated ICU 2.4. This class may be removed or modified. 242 */ 243 CalendarAstronomer::CalendarAstronomer(): 244 fTime(Calendar::getNow()), fLongitude(0.0), fLatitude(0.0), fGmtOffset(0.0), moonPosition(0,0), moonPositionSet(FALSE) { 245 clearCache(); 246 } 247 248 /** 249 * Construct a new <code>CalendarAstronomer</code> object that is initialized to 250 * the specified date and time. 251 * @internal 252 * @deprecated ICU 2.4. This class may be removed or modified. 253 */ 254 CalendarAstronomer::CalendarAstronomer(UDate d): fTime(d), fLongitude(0.0), fLatitude(0.0), fGmtOffset(0.0), moonPosition(0,0), moonPositionSet(FALSE) { 255 clearCache(); 256 } 257 258 /** 259 * Construct a new <code>CalendarAstronomer</code> object with the given 260 * latitude and longitude. The object's time is set to the current 261 * date and time. 262 * <p> 263 * @param longitude The desired longitude, in <em>degrees</em> east of 264 * the Greenwich meridian. 265 * 266 * @param latitude The desired latitude, in <em>degrees</em>. Positive 267 * values signify North, negative South. 268 * 269 * @see java.util.Date#getTime() 270 * @internal 271 * @deprecated ICU 2.4. This class may be removed or modified. 272 */ 273 CalendarAstronomer::CalendarAstronomer(double longitude, double latitude) : 274 fTime(Calendar::getNow()), moonPosition(0,0), moonPositionSet(FALSE) { 275 fLongitude = normPI(longitude * (double)DEG_RAD); 276 fLatitude = normPI(latitude * (double)DEG_RAD); 277 fGmtOffset = (double)(fLongitude * 24. * (double)HOUR_MS / (double)CalendarAstronomer_PI2); 278 clearCache(); 279 } 280 281 CalendarAstronomer::~CalendarAstronomer() 282 { 283 } 284 285 //------------------------------------------------------------------------- 286 // Time and date getters and setters 287 //------------------------------------------------------------------------- 288 289 /** 290 * Set the current date and time of this <code>CalendarAstronomer</code> object. All 291 * astronomical calculations are performed based on this time setting. 292 * 293 * @param aTime the date and time, expressed as the number of milliseconds since 294 * 1/1/1970 0:00 GMT (Gregorian). 295 * 296 * @see #setDate 297 * @see #getTime 298 * @internal 299 * @deprecated ICU 2.4. This class may be removed or modified. 300 */ 301 void CalendarAstronomer::setTime(UDate aTime) { 302 fTime = aTime; 303 U_DEBUG_ASTRO_MSG(("setTime(%.1lf, %sL)\n", aTime, debug_astro_date(aTime+fGmtOffset))); 304 clearCache(); 305 } 306 307 /** 308 * Set the current date and time of this <code>CalendarAstronomer</code> object. All 309 * astronomical calculations are performed based on this time setting. 310 * 311 * @param jdn the desired time, expressed as a "julian day number", 312 * which is the number of elapsed days since 313 * 1/1/4713 BC (Julian), 12:00 GMT. Note that julian day 314 * numbers start at <em>noon</em>. To get the jdn for 315 * the corresponding midnight, subtract 0.5. 316 * 317 * @see #getJulianDay 318 * @see #JULIAN_EPOCH_MS 319 * @internal 320 * @deprecated ICU 2.4. This class may be removed or modified. 321 */ 322 void CalendarAstronomer::setJulianDay(double jdn) { 323 fTime = (double)(jdn * DAY_MS) + JULIAN_EPOCH_MS; 324 clearCache(); 325 julianDay = jdn; 326 } 327 328 /** 329 * Get the current time of this <code>CalendarAstronomer</code> object, 330 * represented as the number of milliseconds since 331 * 1/1/1970 AD 0:00 GMT (Gregorian). 332 * 333 * @see #setTime 334 * @see #getDate 335 * @internal 336 * @deprecated ICU 2.4. This class may be removed or modified. 337 */ 338 UDate CalendarAstronomer::getTime() { 339 return fTime; 340 } 341 342 /** 343 * Get the current time of this <code>CalendarAstronomer</code> object, 344 * expressed as a "julian day number", which is the number of elapsed 345 * days since 1/1/4713 BC (Julian), 12:00 GMT. 346 * 347 * @see #setJulianDay 348 * @see #JULIAN_EPOCH_MS 349 * @internal 350 * @deprecated ICU 2.4. This class may be removed or modified. 351 */ 352 double CalendarAstronomer::getJulianDay() { 353 if (isINVALID(julianDay)) { 354 julianDay = (fTime - (double)JULIAN_EPOCH_MS) / (double)DAY_MS; 355 } 356 return julianDay; 357 } 358 359 /** 360 * Return this object's time expressed in julian centuries: 361 * the number of centuries after 1/1/1900 AD, 12:00 GMT 362 * 363 * @see #getJulianDay 364 * @internal 365 * @deprecated ICU 2.4. This class may be removed or modified. 366 */ 367 double CalendarAstronomer::getJulianCentury() { 368 if (isINVALID(julianCentury)) { 369 julianCentury = (getJulianDay() - 2415020.0) / 36525.0; 370 } 371 return julianCentury; 372 } 373 374 /** 375 * Returns the current Greenwich sidereal time, measured in hours 376 * @internal 377 * @deprecated ICU 2.4. This class may be removed or modified. 378 */ 379 double CalendarAstronomer::getGreenwichSidereal() { 380 if (isINVALID(siderealTime)) { 381 // See page 86 of "Practial Astronomy with your Calculator", 382 // by Peter Duffet-Smith, for details on the algorithm. 383 384 double UT = normalize(fTime/(double)HOUR_MS, 24.); 385 386 siderealTime = normalize(getSiderealOffset() + UT*1.002737909, 24.); 387 } 388 return siderealTime; 389 } 390 391 double CalendarAstronomer::getSiderealOffset() { 392 if (isINVALID(siderealT0)) { 393 double JD = uprv_floor(getJulianDay() - 0.5) + 0.5; 394 double S = JD - 2451545.0; 395 double T = S / 36525.0; 396 siderealT0 = normalize(6.697374558 + 2400.051336*T + 0.000025862*T*T, 24); 397 } 398 return siderealT0; 399 } 400 401 /** 402 * Returns the current local sidereal time, measured in hours 403 * @internal 404 * @deprecated ICU 2.4. This class may be removed or modified. 405 */ 406 double CalendarAstronomer::getLocalSidereal() { 407 return normalize(getGreenwichSidereal() + (fGmtOffset/(double)HOUR_MS), 24.); 408 } 409 410 /** 411 * Converts local sidereal time to Universal Time. 412 * 413 * @param lst The Local Sidereal Time, in hours since sidereal midnight 414 * on this object's current date. 415 * 416 * @return The corresponding Universal Time, in milliseconds since 417 * 1 Jan 1970, GMT. 418 */ 419 double CalendarAstronomer::lstToUT(double lst) { 420 // Convert to local mean time 421 double lt = normalize((lst - getSiderealOffset()) * 0.9972695663, 24); 422 423 // Then find local midnight on this day 424 double base = (DAY_MS * ClockMath::floorDivide(fTime + fGmtOffset,(double)DAY_MS)) - fGmtOffset; 425 426 //out(" lt =" + lt + " hours"); 427 //out(" base=" + new Date(base)); 428 429 return base + (long)(lt * HOUR_MS); 430 } 431 432 433 //------------------------------------------------------------------------- 434 // Coordinate transformations, all based on the current time of this object 435 //------------------------------------------------------------------------- 436 437 /** 438 * Convert from ecliptic to equatorial coordinates. 439 * 440 * @param ecliptic A point in the sky in ecliptic coordinates. 441 * @return The corresponding point in equatorial coordinates. 442 * @internal 443 * @deprecated ICU 2.4. This class may be removed or modified. 444 */ 445 CalendarAstronomer::Equatorial& CalendarAstronomer::eclipticToEquatorial(CalendarAstronomer::Equatorial& result, const CalendarAstronomer::Ecliptic& ecliptic) 446 { 447 return eclipticToEquatorial(result, ecliptic.longitude, ecliptic.latitude); 448 } 449 450 /** 451 * Convert from ecliptic to equatorial coordinates. 452 * 453 * @param eclipLong The ecliptic longitude 454 * @param eclipLat The ecliptic latitude 455 * 456 * @return The corresponding point in equatorial coordinates. 457 * @internal 458 * @deprecated ICU 2.4. This class may be removed or modified. 459 */ 460 CalendarAstronomer::Equatorial& CalendarAstronomer::eclipticToEquatorial(CalendarAstronomer::Equatorial& result, double eclipLong, double eclipLat) 461 { 462 // See page 42 of "Practial Astronomy with your Calculator", 463 // by Peter Duffet-Smith, for details on the algorithm. 464 465 double obliq = eclipticObliquity(); 466 double sinE = ::sin(obliq); 467 double cosE = cos(obliq); 468 469 double sinL = ::sin(eclipLong); 470 double cosL = cos(eclipLong); 471 472 double sinB = ::sin(eclipLat); 473 double cosB = cos(eclipLat); 474 double tanB = tan(eclipLat); 475 476 result.set(atan2(sinL*cosE - tanB*sinE, cosL), 477 asin(sinB*cosE + cosB*sinE*sinL) ); 478 return result; 479 } 480 481 /** 482 * Convert from ecliptic longitude to equatorial coordinates. 483 * 484 * @param eclipLong The ecliptic longitude 485 * 486 * @return The corresponding point in equatorial coordinates. 487 * @internal 488 * @deprecated ICU 2.4. This class may be removed or modified. 489 */ 490 CalendarAstronomer::Equatorial& CalendarAstronomer::eclipticToEquatorial(CalendarAstronomer::Equatorial& result, double eclipLong) 491 { 492 return eclipticToEquatorial(result, eclipLong, 0); // TODO: optimize 493 } 494 495 /** 496 * @internal 497 * @deprecated ICU 2.4. This class may be removed or modified. 498 */ 499 CalendarAstronomer::Horizon& CalendarAstronomer::eclipticToHorizon(CalendarAstronomer::Horizon& result, double eclipLong) 500 { 501 Equatorial equatorial; 502 eclipticToEquatorial(equatorial, eclipLong); 503 504 double H = getLocalSidereal()*CalendarAstronomer::PI/12 - equatorial.ascension; // Hour-angle 505 506 double sinH = ::sin(H); 507 double cosH = cos(H); 508 double sinD = ::sin(equatorial.declination); 509 double cosD = cos(equatorial.declination); 510 double sinL = ::sin(fLatitude); 511 double cosL = cos(fLatitude); 512 513 double altitude = asin(sinD*sinL + cosD*cosL*cosH); 514 double azimuth = atan2(-cosD*cosL*sinH, sinD - sinL * ::sin(altitude)); 515 516 result.set(azimuth, altitude); 517 return result; 518 } 519 520 521 //------------------------------------------------------------------------- 522 // The Sun 523 //------------------------------------------------------------------------- 524 525 // 526 // Parameters of the Sun's orbit as of the epoch Jan 0.0 1990 527 // Angles are in radians (after multiplying by CalendarAstronomer::PI/180) 528 // 529 #define JD_EPOCH 2447891.5 // Julian day of epoch 530 531 #define SUN_ETA_G (279.403303 * CalendarAstronomer::PI/180) // Ecliptic longitude at epoch 532 #define SUN_OMEGA_G (282.768422 * CalendarAstronomer::PI/180) // Ecliptic longitude of perigee 533 #define SUN_E 0.016713 // Eccentricity of orbit 534 //double sunR0 1.495585e8 // Semi-major axis in KM 535 //double sunTheta0 (0.533128 * CalendarAstronomer::PI/180) // Angular diameter at R0 536 537 // The following three methods, which compute the sun parameters 538 // given above for an arbitrary epoch (whatever time the object is 539 // set to), make only a small difference as compared to using the 540 // above constants. E.g., Sunset times might differ by ~12 541 // seconds. Furthermore, the eta-g computation is befuddled by 542 // Duffet-Smith's incorrect coefficients (p.86). I've corrected 543 // the first-order coefficient but the others may be off too - no 544 // way of knowing without consulting another source. 545 546 // /** 547 // * Return the sun's ecliptic longitude at perigee for the current time. 548 // * See Duffett-Smith, p. 86. 549 // * @return radians 550 // */ 551 // private double getSunOmegaG() { 552 // double T = getJulianCentury(); 553 // return (281.2208444 + (1.719175 + 0.000452778*T)*T) * DEG_RAD; 554 // } 555 556 // /** 557 // * Return the sun's ecliptic longitude for the current time. 558 // * See Duffett-Smith, p. 86. 559 // * @return radians 560 // */ 561 // private double getSunEtaG() { 562 // double T = getJulianCentury(); 563 // //return (279.6966778 + (36000.76892 + 0.0003025*T)*T) * DEG_RAD; 564 // // 565 // // The above line is from Duffett-Smith, and yields manifestly wrong 566 // // results. The below constant is derived empirically to match the 567 // // constant he gives for the 1990 EPOCH. 568 // // 569 // return (279.6966778 + (-0.3262541582718024 + 0.0003025*T)*T) * DEG_RAD; 570 // } 571 572 // /** 573 // * Return the sun's eccentricity of orbit for the current time. 574 // * See Duffett-Smith, p. 86. 575 // * @return double 576 // */ 577 // private double getSunE() { 578 // double T = getJulianCentury(); 579 // return 0.01675104 - (0.0000418 + 0.000000126*T)*T; 580 // } 581 582 /** 583 * Find the "true anomaly" (longitude) of an object from 584 * its mean anomaly and the eccentricity of its orbit. This uses 585 * an iterative solution to Kepler's equation. 586 * 587 * @param meanAnomaly The object's longitude calculated as if it were in 588 * a regular, circular orbit, measured in radians 589 * from the point of perigee. 590 * 591 * @param eccentricity The eccentricity of the orbit 592 * 593 * @return The true anomaly (longitude) measured in radians 594 */ 595 static double trueAnomaly(double meanAnomaly, double eccentricity) 596 { 597 // First, solve Kepler's equation iteratively 598 // Duffett-Smith, p.90 599 double delta; 600 double E = meanAnomaly; 601 do { 602 delta = E - eccentricity * ::sin(E) - meanAnomaly; 603 E = E - delta / (1 - eccentricity * ::cos(E)); 604 } 605 while (uprv_fabs(delta) > 1e-5); // epsilon = 1e-5 rad 606 607 return 2.0 * ::atan( ::tan(E/2) * ::sqrt( (1+eccentricity) 608 /(1-eccentricity) ) ); 609 } 610 611 /** 612 * The longitude of the sun at the time specified by this object. 613 * The longitude is measured in radians along the ecliptic 614 * from the "first point of Aries," the point at which the ecliptic 615 * crosses the earth's equatorial plane at the vernal equinox. 616 * <p> 617 * Currently, this method uses an approximation of the two-body Kepler's 618 * equation for the earth and the sun. It does not take into account the 619 * perturbations caused by the other planets, the moon, etc. 620 * @internal 621 * @deprecated ICU 2.4. This class may be removed or modified. 622 */ 623 double CalendarAstronomer::getSunLongitude() 624 { 625 // See page 86 of "Practial Astronomy with your Calculator", 626 // by Peter Duffet-Smith, for details on the algorithm. 627 628 if (isINVALID(sunLongitude)) { 629 getSunLongitude(getJulianDay(), sunLongitude, meanAnomalySun); 630 } 631 return sunLongitude; 632 } 633 634 /** 635 * TODO Make this public when the entire class is package-private. 636 */ 637 /*public*/ void CalendarAstronomer::getSunLongitude(double jDay, double &longitude, double &meanAnomaly) 638 { 639 // See page 86 of "Practial Astronomy with your Calculator", 640 // by Peter Duffet-Smith, for details on the algorithm. 641 642 double day = jDay - JD_EPOCH; // Days since epoch 643 644 // Find the angular distance the sun in a fictitious 645 // circular orbit has travelled since the epoch. 646 double epochAngle = norm2PI(CalendarAstronomer_PI2/TROPICAL_YEAR*day); 647 648 // The epoch wasn't at the sun's perigee; find the angular distance 649 // since perigee, which is called the "mean anomaly" 650 meanAnomaly = norm2PI(epochAngle + SUN_ETA_G - SUN_OMEGA_G); 651 652 // Now find the "true anomaly", e.g. the real solar longitude 653 // by solving Kepler's equation for an elliptical orbit 654 // NOTE: The 3rd ed. of the book lists omega_g and eta_g in different 655 // equations; omega_g is to be correct. 656 longitude = norm2PI(trueAnomaly(meanAnomaly, SUN_E) + SUN_OMEGA_G); 657 } 658 659 /** 660 * The position of the sun at this object's current date and time, 661 * in equatorial coordinates. 662 * @internal 663 * @deprecated ICU 2.4. This class may be removed or modified. 664 */ 665 CalendarAstronomer::Equatorial& CalendarAstronomer::getSunPosition(CalendarAstronomer::Equatorial& result) { 666 return eclipticToEquatorial(result, getSunLongitude(), 0); 667 } 668 669 670 /** 671 * Constant representing the vernal equinox. 672 * For use with {@link #getSunTime getSunTime}. 673 * Note: In this case, "vernal" refers to the northern hemisphere's seasons. 674 * @internal 675 * @deprecated ICU 2.4. This class may be removed or modified. 676 */ 677 /*double CalendarAstronomer::VERNAL_EQUINOX() { 678 return 0; 679 }*/ 680 681 /** 682 * Constant representing the summer solstice. 683 * For use with {@link #getSunTime getSunTime}. 684 * Note: In this case, "summer" refers to the northern hemisphere's seasons. 685 * @internal 686 * @deprecated ICU 2.4. This class may be removed or modified. 687 */ 688 double CalendarAstronomer::SUMMER_SOLSTICE() { 689 return (CalendarAstronomer::PI/2); 690 } 691 692 /** 693 * Constant representing the autumnal equinox. 694 * For use with {@link #getSunTime getSunTime}. 695 * Note: In this case, "autumn" refers to the northern hemisphere's seasons. 696 * @internal 697 * @deprecated ICU 2.4. This class may be removed or modified. 698 */ 699 /*double CalendarAstronomer::AUTUMN_EQUINOX() { 700 return (CalendarAstronomer::PI); 701 }*/ 702 703 /** 704 * Constant representing the winter solstice. 705 * For use with {@link #getSunTime getSunTime}. 706 * Note: In this case, "winter" refers to the northern hemisphere's seasons. 707 * @internal 708 * @deprecated ICU 2.4. This class may be removed or modified. 709 */ 710 double CalendarAstronomer::WINTER_SOLSTICE() { 711 return ((CalendarAstronomer::PI*3)/2); 712 } 713 714 CalendarAstronomer::AngleFunc::~AngleFunc() {} 715 716 /** 717 * Find the next time at which the sun's ecliptic longitude will have 718 * the desired value. 719 * @internal 720 * @deprecated ICU 2.4. This class may be removed or modified. 721 */ 722 class SunTimeAngleFunc : public CalendarAstronomer::AngleFunc { 723 public: 724 virtual double eval(CalendarAstronomer& a) { return a.getSunLongitude(); } 725 }; 726 727 UDate CalendarAstronomer::getSunTime(double desired, UBool next) 728 { 729 SunTimeAngleFunc func; 730 return timeOfAngle( func, 731 desired, 732 TROPICAL_YEAR, 733 MINUTE_MS, 734 next); 735 } 736 737 CalendarAstronomer::CoordFunc::~CoordFunc() {} 738 739 class RiseSetCoordFunc : public CalendarAstronomer::CoordFunc { 740 public: 741 virtual void eval(CalendarAstronomer::Equatorial& result, CalendarAstronomer&a) { a.getSunPosition(result); } 742 }; 743 744 UDate CalendarAstronomer::getSunRiseSet(UBool rise) 745 { 746 UDate t0 = fTime; 747 748 // Make a rough guess: 6am or 6pm local time on the current day 749 double noon = ClockMath::floorDivide(fTime + fGmtOffset, (double)DAY_MS)*DAY_MS - fGmtOffset + (12*HOUR_MS); 750 751 U_DEBUG_ASTRO_MSG(("Noon=%.2lf, %sL, gmtoff %.2lf\n", noon, debug_astro_date(noon+fGmtOffset), fGmtOffset)); 752 setTime(noon + ((rise ? -6 : 6) * HOUR_MS)); 753 U_DEBUG_ASTRO_MSG(("added %.2lf ms as a guess,\n", ((rise ? -6. : 6.) * HOUR_MS))); 754 755 RiseSetCoordFunc func; 756 double t = riseOrSet(func, 757 rise, 758 .533 * DEG_RAD, // Angular Diameter 759 34. /60.0 * DEG_RAD, // Refraction correction 760 MINUTE_MS / 12.); // Desired accuracy 761 762 setTime(t0); 763 return t; 764 } 765 766 // Commented out - currently unused. ICU 2.6, Alan 767 // //------------------------------------------------------------------------- 768 // // Alternate Sun Rise/Set 769 // // See Duffett-Smith p.93 770 // //------------------------------------------------------------------------- 771 // 772 // // This yields worse results (as compared to USNO data) than getSunRiseSet(). 773 // /** 774 // * TODO Make this when the entire class is package-private. 775 // */ 776 // /*public*/ long getSunRiseSet2(boolean rise) { 777 // // 1. Calculate coordinates of the sun's center for midnight 778 // double jd = uprv_floor(getJulianDay() - 0.5) + 0.5; 779 // double[] sl = getSunLongitude(jd);// double lambda1 = sl[0]; 780 // Equatorial pos1 = eclipticToEquatorial(lambda1, 0); 781 // 782 // // 2. Add ... to lambda to get position 24 hours later 783 // double lambda2 = lambda1 + 0.985647*DEG_RAD; 784 // Equatorial pos2 = eclipticToEquatorial(lambda2, 0); 785 // 786 // // 3. Calculate LSTs of rising and setting for these two positions 787 // double tanL = ::tan(fLatitude); 788 // double H = ::acos(-tanL * ::tan(pos1.declination)); 789 // double lst1r = (CalendarAstronomer_PI2 + pos1.ascension - H) * 24 / CalendarAstronomer_PI2; 790 // double lst1s = (pos1.ascension + H) * 24 / CalendarAstronomer_PI2; 791 // H = ::acos(-tanL * ::tan(pos2.declination)); 792 // double lst2r = (CalendarAstronomer_PI2-H + pos2.ascension ) * 24 / CalendarAstronomer_PI2; 793 // double lst2s = (H + pos2.ascension ) * 24 / CalendarAstronomer_PI2; 794 // if (lst1r > 24) lst1r -= 24; 795 // if (lst1s > 24) lst1s -= 24; 796 // if (lst2r > 24) lst2r -= 24; 797 // if (lst2s > 24) lst2s -= 24; 798 // 799 // // 4. Convert LSTs to GSTs. If GST1 > GST2, add 24 to GST2. 800 // double gst1r = lstToGst(lst1r); 801 // double gst1s = lstToGst(lst1s); 802 // double gst2r = lstToGst(lst2r); 803 // double gst2s = lstToGst(lst2s); 804 // if (gst1r > gst2r) gst2r += 24; 805 // if (gst1s > gst2s) gst2s += 24; 806 // 807 // // 5. Calculate GST at 0h UT of this date 808 // double t00 = utToGst(0); 809 // 810 // // 6. Calculate GST at 0h on the observer's longitude 811 // double offset = ::round(fLongitude*12/PI); // p.95 step 6; he _rounds_ to nearest 15 deg. 812 // double t00p = t00 - offset*1.002737909; 813 // if (t00p < 0) t00p += 24; // do NOT normalize 814 // 815 // // 7. Adjust 816 // if (gst1r < t00p) { 817 // gst1r += 24; 818 // gst2r += 24; 819 // } 820 // if (gst1s < t00p) { 821 // gst1s += 24; 822 // gst2s += 24; 823 // } 824 // 825 // // 8. 826 // double gstr = (24.07*gst1r-t00*(gst2r-gst1r))/(24.07+gst1r-gst2r); 827 // double gsts = (24.07*gst1s-t00*(gst2s-gst1s))/(24.07+gst1s-gst2s); 828 // 829 // // 9. Correct for parallax, refraction, and sun's diameter 830 // double dec = (pos1.declination + pos2.declination) / 2; 831 // double psi = ::acos(sin(fLatitude) / cos(dec)); 832 // double x = 0.830725 * DEG_RAD; // parallax+refraction+diameter 833 // double y = ::asin(sin(x) / ::sin(psi)) * RAD_DEG; 834 // double delta_t = 240 * y / cos(dec) / 3600; // hours 835 // 836 // // 10. Add correction to GSTs, subtract from GSTr 837 // gstr -= delta_t; 838 // gsts += delta_t; 839 // 840 // // 11. Convert GST to UT and then to local civil time 841 // double ut = gstToUt(rise ? gstr : gsts); 842 // //System.out.println((rise?"rise=":"set=") + ut + ", delta_t=" + delta_t); 843 // long midnight = DAY_MS * (time / DAY_MS); // Find UT midnight on this day 844 // return midnight + (long) (ut * 3600000); 845 // } 846 847 // Commented out - currently unused. ICU 2.6, Alan 848 // /** 849 // * Convert local sidereal time to Greenwich sidereal time. 850 // * Section 15. Duffett-Smith p.21 851 // * @param lst in hours (0..24) 852 // * @return GST in hours (0..24) 853 // */ 854 // double lstToGst(double lst) { 855 // double delta = fLongitude * 24 / CalendarAstronomer_PI2; 856 // return normalize(lst - delta, 24); 857 // } 858 859 // Commented out - currently unused. ICU 2.6, Alan 860 // /** 861 // * Convert UT to GST on this date. 862 // * Section 12. Duffett-Smith p.17 863 // * @param ut in hours 864 // * @return GST in hours 865 // */ 866 // double utToGst(double ut) { 867 // return normalize(getT0() + ut*1.002737909, 24); 868 // } 869 870 // Commented out - currently unused. ICU 2.6, Alan 871 // /** 872 // * Convert GST to UT on this date. 873 // * Section 13. Duffett-Smith p.18 874 // * @param gst in hours 875 // * @return UT in hours 876 // */ 877 // double gstToUt(double gst) { 878 // return normalize(gst - getT0(), 24) * 0.9972695663; 879 // } 880 881 // Commented out - currently unused. ICU 2.6, Alan 882 // double getT0() { 883 // // Common computation for UT <=> GST 884 // 885 // // Find JD for 0h UT 886 // double jd = uprv_floor(getJulianDay() - 0.5) + 0.5; 887 // 888 // double s = jd - 2451545.0; 889 // double t = s / 36525.0; 890 // double t0 = 6.697374558 + (2400.051336 + 0.000025862*t)*t; 891 // return t0; 892 // } 893 894 // Commented out - currently unused. ICU 2.6, Alan 895 // //------------------------------------------------------------------------- 896 // // Alternate Sun Rise/Set 897 // // See sci.astro FAQ 898 // // http://www.faqs.org/faqs/astronomy/faq/part3/section-5.html 899 // //------------------------------------------------------------------------- 900 // 901 // // Note: This method appears to produce inferior accuracy as 902 // // compared to getSunRiseSet(). 903 // 904 // /** 905 // * TODO Make this when the entire class is package-private. 906 // */ 907 // /*public*/ long getSunRiseSet3(boolean rise) { 908 // 909 // // Compute day number for 0.0 Jan 2000 epoch 910 // double d = (double)(time - EPOCH_2000_MS) / DAY_MS; 911 // 912 // // Now compute the Local Sidereal Time, LST: 913 // // 914 // double LST = 98.9818 + 0.985647352 * d + /*UT*15 + long*/ 915 // fLongitude*RAD_DEG; 916 // // 917 // // (east long. positive). Note that LST is here expressed in degrees, 918 // // where 15 degrees corresponds to one hour. Since LST really is an angle, 919 // // it's convenient to use one unit---degrees---throughout. 920 // 921 // // COMPUTING THE SUN'S POSITION 922 // // ---------------------------- 923 // // 924 // // To be able to compute the Sun's rise/set times, you need to be able to 925 // // compute the Sun's position at any time. First compute the "day 926 // // number" d as outlined above, for the desired moment. Next compute: 927 // // 928 // double oblecl = 23.4393 - 3.563E-7 * d; 929 // // 930 // double w = 282.9404 + 4.70935E-5 * d; 931 // double M = 356.0470 + 0.9856002585 * d; 932 // double e = 0.016709 - 1.151E-9 * d; 933 // // 934 // // This is the obliquity of the ecliptic, plus some of the elements of 935 // // the Sun's apparent orbit (i.e., really the Earth's orbit): w = 936 // // argument of perihelion, M = mean anomaly, e = eccentricity. 937 // // Semi-major axis is here assumed to be exactly 1.0 (while not strictly 938 // // true, this is still an accurate approximation). Next compute E, the 939 // // eccentric anomaly: 940 // // 941 // double E = M + e*(180/PI) * ::sin(M*DEG_RAD) * ( 1.0 + e*cos(M*DEG_RAD) ); 942 // // 943 // // where E and M are in degrees. This is it---no further iterations are 944 // // needed because we know e has a sufficiently small value. Next compute 945 // // the true anomaly, v, and the distance, r: 946 // // 947 // /* r * cos(v) = */ double A = cos(E*DEG_RAD) - e; 948 // /* r * ::sin(v) = */ double B = ::sqrt(1 - e*e) * ::sin(E*DEG_RAD); 949 // // 950 // // and 951 // // 952 // // r = sqrt( A*A + B*B ) 953 // double v = ::atan2( B, A )*RAD_DEG; 954 // // 955 // // The Sun's true longitude, slon, can now be computed: 956 // // 957 // double slon = v + w; 958 // // 959 // // Since the Sun is always at the ecliptic (or at least very very close to 960 // // it), we can use simplified formulae to convert slon (the Sun's ecliptic 961 // // longitude) to sRA and sDec (the Sun's RA and Dec): 962 // // 963 // // ::sin(slon) * cos(oblecl) 964 // // tan(sRA) = ------------------------- 965 // // cos(slon) 966 // // 967 // // ::sin(sDec) = ::sin(oblecl) * ::sin(slon) 968 // // 969 // // As was the case when computing az, the Azimuth, if possible use an 970 // // atan2() function to compute sRA. 971 // 972 // double sRA = ::atan2(sin(slon*DEG_RAD) * cos(oblecl*DEG_RAD), cos(slon*DEG_RAD))*RAD_DEG; 973 // 974 // double sin_sDec = ::sin(oblecl*DEG_RAD) * ::sin(slon*DEG_RAD); 975 // double sDec = ::asin(sin_sDec)*RAD_DEG; 976 // 977 // // COMPUTING RISE AND SET TIMES 978 // // ---------------------------- 979 // // 980 // // To compute when an object rises or sets, you must compute when it 981 // // passes the meridian and the HA of rise/set. Then the rise time is 982 // // the meridian time minus HA for rise/set, and the set time is the 983 // // meridian time plus the HA for rise/set. 984 // // 985 // // To find the meridian time, compute the Local Sidereal Time at 0h local 986 // // time (or 0h UT if you prefer to work in UT) as outlined above---name 987 // // that quantity LST0. The Meridian Time, MT, will now be: 988 // // 989 // // MT = RA - LST0 990 // double MT = normalize(sRA - LST, 360); 991 // // 992 // // where "RA" is the object's Right Ascension (in degrees!). If negative, 993 // // add 360 deg to MT. If the object is the Sun, leave the time as it is, 994 // // but if it's stellar, multiply MT by 365.2422/366.2422, to convert from 995 // // sidereal to solar time. Now, compute HA for rise/set, name that 996 // // quantity HA0: 997 // // 998 // // ::sin(h0) - ::sin(lat) * ::sin(Dec) 999 // // cos(HA0) = --------------------------------- 1000 // // cos(lat) * cos(Dec) 1001 // // 1002 // // where h0 is the altitude selected to represent rise/set. For a purely 1003 // // mathematical horizon, set h0 = 0 and simplify to: 1004 // // 1005 // // cos(HA0) = - tan(lat) * tan(Dec) 1006 // // 1007 // // If you want to account for refraction on the atmosphere, set h0 = -35/60 1008 // // degrees (-35 arc minutes), and if you want to compute the rise/set times 1009 // // for the Sun's upper limb, set h0 = -50/60 (-50 arc minutes). 1010 // // 1011 // double h0 = -50/60 * DEG_RAD; 1012 // 1013 // double HA0 = ::acos( 1014 // (sin(h0) - ::sin(fLatitude) * sin_sDec) / 1015 // (cos(fLatitude) * cos(sDec*DEG_RAD)))*RAD_DEG; 1016 // 1017 // // When HA0 has been computed, leave it as it is for the Sun but multiply 1018 // // by 365.2422/366.2422 for stellar objects, to convert from sidereal to 1019 // // solar time. Finally compute: 1020 // // 1021 // // Rise time = MT - HA0 1022 // // Set time = MT + HA0 1023 // // 1024 // // convert the times from degrees to hours by dividing by 15. 1025 // // 1026 // // If you'd like to check that your calculations are accurate or just 1027 // // need a quick result, check the USNO's Sun or Moon Rise/Set Table, 1028 // // <URL:http://aa.usno.navy.mil/AA/data/docs/RS_OneYear.html>. 1029 // 1030 // double result = MT + (rise ? -HA0 : HA0); // in degrees 1031 // 1032 // // Find UT midnight on this day 1033 // long midnight = DAY_MS * (time / DAY_MS); 1034 // 1035 // return midnight + (long) (result * 3600000 / 15); 1036 // } 1037 1038 //------------------------------------------------------------------------- 1039 // The Moon 1040 //------------------------------------------------------------------------- 1041 1042 #define moonL0 (318.351648 * CalendarAstronomer::PI/180 ) // Mean long. at epoch 1043 #define moonP0 ( 36.340410 * CalendarAstronomer::PI/180 ) // Mean long. of perigee 1044 #define moonN0 ( 318.510107 * CalendarAstronomer::PI/180 ) // Mean long. of node 1045 #define moonI ( 5.145366 * CalendarAstronomer::PI/180 ) // Inclination of orbit 1046 #define moonE ( 0.054900 ) // Eccentricity of orbit 1047 1048 // These aren't used right now 1049 #define moonA ( 3.84401e5 ) // semi-major axis (km) 1050 #define moonT0 ( 0.5181 * CalendarAstronomer::PI/180 ) // Angular size at distance A 1051 #define moonPi ( 0.9507 * CalendarAstronomer::PI/180 ) // Parallax at distance A 1052 1053 /** 1054 * The position of the moon at the time set on this 1055 * object, in equatorial coordinates. 1056 * @internal 1057 * @deprecated ICU 2.4. This class may be removed or modified. 1058 */ 1059 const CalendarAstronomer::Equatorial& CalendarAstronomer::getMoonPosition() 1060 { 1061 // 1062 // See page 142 of "Practial Astronomy with your Calculator", 1063 // by Peter Duffet-Smith, for details on the algorithm. 1064 // 1065 if (moonPositionSet == FALSE) { 1066 // Calculate the solar longitude. Has the side effect of 1067 // filling in "meanAnomalySun" as well. 1068 getSunLongitude(); 1069 1070 // 1071 // Find the # of days since the epoch of our orbital parameters. 1072 // TODO: Convert the time of day portion into ephemeris time 1073 // 1074 double day = getJulianDay() - JD_EPOCH; // Days since epoch 1075 1076 // Calculate the mean longitude and anomaly of the moon, based on 1077 // a circular orbit. Similar to the corresponding solar calculation. 1078 double meanLongitude = norm2PI(13.1763966*PI/180*day + moonL0); 1079 meanAnomalyMoon = norm2PI(meanLongitude - 0.1114041*PI/180 * day - moonP0); 1080 1081 // 1082 // Calculate the following corrections: 1083 // Evection: the sun's gravity affects the moon's eccentricity 1084 // Annual Eqn: variation in the effect due to earth-sun distance 1085 // A3: correction factor (for ???) 1086 // 1087 double evection = 1.2739*PI/180 * ::sin(2 * (meanLongitude - sunLongitude) 1088 - meanAnomalyMoon); 1089 double annual = 0.1858*PI/180 * ::sin(meanAnomalySun); 1090 double a3 = 0.3700*PI/180 * ::sin(meanAnomalySun); 1091 1092 meanAnomalyMoon += evection - annual - a3; 1093 1094 // 1095 // More correction factors: 1096 // center equation of the center correction 1097 // a4 yet another error correction (???) 1098 // 1099 // TODO: Skip the equation of the center correction and solve Kepler's eqn? 1100 // 1101 double center = 6.2886*PI/180 * ::sin(meanAnomalyMoon); 1102 double a4 = 0.2140*PI/180 * ::sin(2 * meanAnomalyMoon); 1103 1104 // Now find the moon's corrected longitude 1105 moonLongitude = meanLongitude + evection + center - annual + a4; 1106 1107 // 1108 // And finally, find the variation, caused by the fact that the sun's 1109 // gravitational pull on the moon varies depending on which side of 1110 // the earth the moon is on 1111 // 1112 double variation = 0.6583*CalendarAstronomer::PI/180 * ::sin(2*(moonLongitude - sunLongitude)); 1113 1114 moonLongitude += variation; 1115 1116 // 1117 // What we've calculated so far is the moon's longitude in the plane 1118 // of its own orbit. Now map to the ecliptic to get the latitude 1119 // and longitude. First we need to find the longitude of the ascending 1120 // node, the position on the ecliptic where it is crossed by the moon's 1121 // orbit as it crosses from the southern to the northern hemisphere. 1122 // 1123 double nodeLongitude = norm2PI(moonN0 - 0.0529539*PI/180 * day); 1124 1125 nodeLongitude -= 0.16*PI/180 * ::sin(meanAnomalySun); 1126 1127 double y = ::sin(moonLongitude - nodeLongitude); 1128 double x = cos(moonLongitude - nodeLongitude); 1129 1130 moonEclipLong = ::atan2(y*cos(moonI), x) + nodeLongitude; 1131 double moonEclipLat = ::asin(y * ::sin(moonI)); 1132 1133 eclipticToEquatorial(moonPosition, moonEclipLong, moonEclipLat); 1134 moonPositionSet = TRUE; 1135 } 1136 return moonPosition; 1137 } 1138 1139 /** 1140 * The "age" of the moon at the time specified in this object. 1141 * This is really the angle between the 1142 * current ecliptic longitudes of the sun and the moon, 1143 * measured in radians. 1144 * 1145 * @see #getMoonPhase 1146 * @internal 1147 * @deprecated ICU 2.4. This class may be removed or modified. 1148 */ 1149 double CalendarAstronomer::getMoonAge() { 1150 // See page 147 of "Practial Astronomy with your Calculator", 1151 // by Peter Duffet-Smith, for details on the algorithm. 1152 // 1153 // Force the moon's position to be calculated. We're going to use 1154 // some the intermediate results cached during that calculation. 1155 // 1156 getMoonPosition(); 1157 1158 return norm2PI(moonEclipLong - sunLongitude); 1159 } 1160 1161 /** 1162 * Calculate the phase of the moon at the time set in this object. 1163 * The returned phase is a <code>double</code> in the range 1164 * <code>0 <= phase < 1</code>, interpreted as follows: 1165 * <ul> 1166 * <li>0.00: New moon 1167 * <li>0.25: First quarter 1168 * <li>0.50: Full moon 1169 * <li>0.75: Last quarter 1170 * </ul> 1171 * 1172 * @see #getMoonAge 1173 * @internal 1174 * @deprecated ICU 2.4. This class may be removed or modified. 1175 */ 1176 double CalendarAstronomer::getMoonPhase() { 1177 // See page 147 of "Practial Astronomy with your Calculator", 1178 // by Peter Duffet-Smith, for details on the algorithm. 1179 return 0.5 * (1 - cos(getMoonAge())); 1180 } 1181 1182 /** 1183 * Constant representing a new moon. 1184 * For use with {@link #getMoonTime getMoonTime} 1185 * @internal 1186 * @deprecated ICU 2.4. This class may be removed or modified. 1187 */ 1188 const CalendarAstronomer::MoonAge CalendarAstronomer::NEW_MOON() { 1189 return CalendarAstronomer::MoonAge(0); 1190 } 1191 1192 /** 1193 * Constant representing the moon's first quarter. 1194 * For use with {@link #getMoonTime getMoonTime} 1195 * @internal 1196 * @deprecated ICU 2.4. This class may be removed or modified. 1197 */ 1198 /*const CalendarAstronomer::MoonAge CalendarAstronomer::FIRST_QUARTER() { 1199 return CalendarAstronomer::MoonAge(CalendarAstronomer::PI/2); 1200 }*/ 1201 1202 /** 1203 * Constant representing a full moon. 1204 * For use with {@link #getMoonTime getMoonTime} 1205 * @internal 1206 * @deprecated ICU 2.4. This class may be removed or modified. 1207 */ 1208 const CalendarAstronomer::MoonAge CalendarAstronomer::FULL_MOON() { 1209 return CalendarAstronomer::MoonAge(CalendarAstronomer::PI); 1210 } 1211 /** 1212 * Constant representing the moon's last quarter. 1213 * For use with {@link #getMoonTime getMoonTime} 1214 * @internal 1215 * @deprecated ICU 2.4. This class may be removed or modified. 1216 */ 1217 1218 class MoonTimeAngleFunc : public CalendarAstronomer::AngleFunc { 1219 public: 1220 virtual double eval(CalendarAstronomer&a) { return a.getMoonAge(); } 1221 }; 1222 1223 /*const CalendarAstronomer::MoonAge CalendarAstronomer::LAST_QUARTER() { 1224 return CalendarAstronomer::MoonAge((CalendarAstronomer::PI*3)/2); 1225 }*/ 1226 1227 /** 1228 * Find the next or previous time at which the Moon's ecliptic 1229 * longitude will have the desired value. 1230 * <p> 1231 * @param desired The desired longitude. 1232 * @param next <tt>true</tt> if the next occurrance of the phase 1233 * is desired, <tt>false</tt> for the previous occurrance. 1234 * @internal 1235 * @deprecated ICU 2.4. This class may be removed or modified. 1236 */ 1237 UDate CalendarAstronomer::getMoonTime(double desired, UBool next) 1238 { 1239 MoonTimeAngleFunc func; 1240 return timeOfAngle( func, 1241 desired, 1242 SYNODIC_MONTH, 1243 MINUTE_MS, 1244 next); 1245 } 1246 1247 /** 1248 * Find the next or previous time at which the moon will be in the 1249 * desired phase. 1250 * <p> 1251 * @param desired The desired phase of the moon. 1252 * @param next <tt>true</tt> if the next occurrance of the phase 1253 * is desired, <tt>false</tt> for the previous occurrance. 1254 * @internal 1255 * @deprecated ICU 2.4. This class may be removed or modified. 1256 */ 1257 UDate CalendarAstronomer::getMoonTime(const CalendarAstronomer::MoonAge& desired, UBool next) { 1258 return getMoonTime(desired.value, next); 1259 } 1260 1261 class MoonRiseSetCoordFunc : public CalendarAstronomer::CoordFunc { 1262 public: 1263 virtual void eval(CalendarAstronomer::Equatorial& result, CalendarAstronomer&a) { result = a.getMoonPosition(); } 1264 }; 1265 1266 /** 1267 * Returns the time (GMT) of sunrise or sunset on the local date to which 1268 * this calendar is currently set. 1269 * @internal 1270 * @deprecated ICU 2.4. This class may be removed or modified. 1271 */ 1272 UDate CalendarAstronomer::getMoonRiseSet(UBool rise) 1273 { 1274 MoonRiseSetCoordFunc func; 1275 return riseOrSet(func, 1276 rise, 1277 .533 * DEG_RAD, // Angular Diameter 1278 34 /60.0 * DEG_RAD, // Refraction correction 1279 MINUTE_MS); // Desired accuracy 1280 } 1281 1282 //------------------------------------------------------------------------- 1283 // Interpolation methods for finding the time at which a given event occurs 1284 //------------------------------------------------------------------------- 1285 1286 UDate CalendarAstronomer::timeOfAngle(AngleFunc& func, double desired, 1287 double periodDays, double epsilon, UBool next) 1288 { 1289 // Find the value of the function at the current time 1290 double lastAngle = func.eval(*this); 1291 1292 // Find out how far we are from the desired angle 1293 double deltaAngle = norm2PI(desired - lastAngle) ; 1294 1295 // Using the average period, estimate the next (or previous) time at 1296 // which the desired angle occurs. 1297 double deltaT = (deltaAngle + (next ? 0.0 : - CalendarAstronomer_PI2 )) * (periodDays*DAY_MS) / CalendarAstronomer_PI2; 1298 1299 double lastDeltaT = deltaT; // Liu 1300 UDate startTime = fTime; // Liu 1301 1302 setTime(fTime + uprv_ceil(deltaT)); 1303 1304 // Now iterate until we get the error below epsilon. Throughout 1305 // this loop we use normPI to get values in the range -Pi to Pi, 1306 // since we're using them as correction factors rather than absolute angles. 1307 do { 1308 // Evaluate the function at the time we've estimated 1309 double angle = func.eval(*this); 1310 1311 // Find the # of milliseconds per radian at this point on the curve 1312 double factor = uprv_fabs(deltaT / normPI(angle-lastAngle)); 1313 1314 // Correct the time estimate based on how far off the angle is 1315 deltaT = normPI(desired - angle) * factor; 1316 1317 // HACK: 1318 // 1319 // If abs(deltaT) begins to diverge we need to quit this loop. 1320 // This only appears to happen when attempting to locate, for 1321 // example, a new moon on the day of the new moon. E.g.: 1322 // 1323 // This result is correct: 1324 // newMoon(7508(Mon Jul 23 00:00:00 CST 1990,false))= 1325 // Sun Jul 22 10:57:41 CST 1990 1326 // 1327 // But attempting to make the same call a day earlier causes deltaT 1328 // to diverge: 1329 // CalendarAstronomer.timeOfAngle() diverging: 1.348508727575625E9 -> 1330 // 1.3649828540224032E9 1331 // newMoon(7507(Sun Jul 22 00:00:00 CST 1990,false))= 1332 // Sun Jul 08 13:56:15 CST 1990 1333 // 1334 // As a temporary solution, we catch this specific condition and 1335 // adjust our start time by one eighth period days (either forward 1336 // or backward) and try again. 1337 // Liu 11/9/00 1338 if (uprv_fabs(deltaT) > uprv_fabs(lastDeltaT)) { 1339 double delta = uprv_ceil (periodDays * DAY_MS / 8.0); 1340 setTime(startTime + (next ? delta : -delta)); 1341 return timeOfAngle(func, desired, periodDays, epsilon, next); 1342 } 1343 1344 lastDeltaT = deltaT; 1345 lastAngle = angle; 1346 1347 setTime(fTime + uprv_ceil(deltaT)); 1348 } 1349 while (uprv_fabs(deltaT) > epsilon); 1350 1351 return fTime; 1352 } 1353 1354 UDate CalendarAstronomer::riseOrSet(CoordFunc& func, UBool rise, 1355 double diameter, double refraction, 1356 double epsilon) 1357 { 1358 Equatorial pos; 1359 double tanL = ::tan(fLatitude); 1360 double deltaT = 0; 1361 int32_t count = 0; 1362 1363 // 1364 // Calculate the object's position at the current time, then use that 1365 // position to calculate the time of rising or setting. The position 1366 // will be different at that time, so iterate until the error is allowable. 1367 // 1368 U_DEBUG_ASTRO_MSG(("setup rise=%s, dia=%.3lf, ref=%.3lf, eps=%.3lf\n", 1369 rise?"T":"F", diameter, refraction, epsilon)); 1370 do { 1371 // See "Practical Astronomy With Your Calculator, section 33. 1372 func.eval(pos, *this); 1373 double angle = ::acos(-tanL * ::tan(pos.declination)); 1374 double lst = ((rise ? CalendarAstronomer_PI2-angle : angle) + pos.ascension ) * 24 / CalendarAstronomer_PI2; 1375 1376 // Convert from LST to Universal Time. 1377 UDate newTime = lstToUT( lst ); 1378 1379 deltaT = newTime - fTime; 1380 setTime(newTime); 1381 U_DEBUG_ASTRO_MSG(("%d] dT=%.3lf, angle=%.3lf, lst=%.3lf, A=%.3lf/D=%.3lf\n", 1382 count, deltaT, angle, lst, pos.ascension, pos.declination)); 1383 } 1384 while (++ count < 5 && uprv_fabs(deltaT) > epsilon); 1385 1386 // Calculate the correction due to refraction and the object's angular diameter 1387 double cosD = ::cos(pos.declination); 1388 double psi = ::acos(sin(fLatitude) / cosD); 1389 double x = diameter / 2 + refraction; 1390 double y = ::asin(sin(x) / ::sin(psi)); 1391 long delta = (long)((240 * y * RAD_DEG / cosD)*SECOND_MS); 1392 1393 return fTime + (rise ? -delta : delta); 1394 } 1395 /** 1396 * Return the obliquity of the ecliptic (the angle between the ecliptic 1397 * and the earth's equator) at the current time. This varies due to 1398 * the precession of the earth's axis. 1399 * 1400 * @return the obliquity of the ecliptic relative to the equator, 1401 * measured in radians. 1402 */ 1403 double CalendarAstronomer::eclipticObliquity() { 1404 if (isINVALID(eclipObliquity)) { 1405 const double epoch = 2451545.0; // 2000 AD, January 1.5 1406 1407 double T = (getJulianDay() - epoch) / 36525; 1408 1409 eclipObliquity = 23.439292 1410 - 46.815/3600 * T 1411 - 0.0006/3600 * T*T 1412 + 0.00181/3600 * T*T*T; 1413 1414 eclipObliquity *= DEG_RAD; 1415 } 1416 return eclipObliquity; 1417 } 1418 1419 1420 //------------------------------------------------------------------------- 1421 // Private data 1422 //------------------------------------------------------------------------- 1423 void CalendarAstronomer::clearCache() { 1424 const double INVALID = uprv_getNaN(); 1425 1426 julianDay = INVALID; 1427 julianCentury = INVALID; 1428 sunLongitude = INVALID; 1429 meanAnomalySun = INVALID; 1430 moonLongitude = INVALID; 1431 moonEclipLong = INVALID; 1432 meanAnomalyMoon = INVALID; 1433 eclipObliquity = INVALID; 1434 siderealTime = INVALID; 1435 siderealT0 = INVALID; 1436 moonPositionSet = FALSE; 1437 } 1438 1439 //private static void out(String s) { 1440 // System.out.println(s); 1441 //} 1442 1443 //private static String deg(double rad) { 1444 // return Double.toString(rad * RAD_DEG); 1445 //} 1446 1447 //private static String hours(long ms) { 1448 // return Double.toString((double)ms / HOUR_MS) + " hours"; 1449 //} 1450 1451 /** 1452 * @internal 1453 * @deprecated ICU 2.4. This class may be removed or modified. 1454 */ 1455 /*UDate CalendarAstronomer::local(UDate localMillis) { 1456 // TODO - srl ? 1457 TimeZone *tz = TimeZone::createDefault(); 1458 int32_t rawOffset; 1459 int32_t dstOffset; 1460 UErrorCode status = U_ZERO_ERROR; 1461 tz->getOffset(localMillis, TRUE, rawOffset, dstOffset, status); 1462 delete tz; 1463 return localMillis - rawOffset; 1464 }*/ 1465 1466 // Debugging functions 1467 UnicodeString CalendarAstronomer::Ecliptic::toString() const 1468 { 1469 #ifdef U_DEBUG_ASTRO 1470 char tmp[800]; 1471 sprintf(tmp, "[%.5f,%.5f]", longitude*RAD_DEG, latitude*RAD_DEG); 1472 return UnicodeString(tmp, ""); 1473 #else 1474 return UnicodeString(); 1475 #endif 1476 } 1477 1478 UnicodeString CalendarAstronomer::Equatorial::toString() const 1479 { 1480 #ifdef U_DEBUG_ASTRO 1481 char tmp[400]; 1482 sprintf(tmp, "%f,%f", 1483 (ascension*RAD_DEG), (declination*RAD_DEG)); 1484 return UnicodeString(tmp, ""); 1485 #else 1486 return UnicodeString(); 1487 #endif 1488 } 1489 1490 UnicodeString CalendarAstronomer::Horizon::toString() const 1491 { 1492 #ifdef U_DEBUG_ASTRO 1493 char tmp[800]; 1494 sprintf(tmp, "[%.5f,%.5f]", altitude*RAD_DEG, azimuth*RAD_DEG); 1495 return UnicodeString(tmp, ""); 1496 #else 1497 return UnicodeString(); 1498 #endif 1499 } 1500 1501 1502 // static private String radToHms(double angle) { 1503 // int hrs = (int) (angle*RAD_HOUR); 1504 // int min = (int)((angle*RAD_HOUR - hrs) * 60); 1505 // int sec = (int)((angle*RAD_HOUR - hrs - min/60.0) * 3600); 1506 1507 // return Integer.toString(hrs) + "h" + min + "m" + sec + "s"; 1508 // } 1509 1510 // static private String radToDms(double angle) { 1511 // int deg = (int) (angle*RAD_DEG); 1512 // int min = (int)((angle*RAD_DEG - deg) * 60); 1513 // int sec = (int)((angle*RAD_DEG - deg - min/60.0) * 3600); 1514 1515 // return Integer.toString(deg) + "\u00b0" + min + "'" + sec + "\""; 1516 // } 1517 1518 // =============== Calendar Cache ================ 1519 1520 void CalendarCache::createCache(CalendarCache** cache, UErrorCode& status) { 1521 ucln_i18n_registerCleanup(UCLN_I18N_ASTRO_CALENDAR, calendar_astro_cleanup); 1522 if(cache == NULL) { 1523 status = U_MEMORY_ALLOCATION_ERROR; 1524 } else { 1525 *cache = new CalendarCache(32, status); 1526 if(U_FAILURE(status)) { 1527 delete *cache; 1528 *cache = NULL; 1529 } 1530 } 1531 } 1532 1533 int32_t CalendarCache::get(CalendarCache** cache, int32_t key, UErrorCode &status) { 1534 int32_t res; 1535 1536 if(U_FAILURE(status)) { 1537 return 0; 1538 } 1539 umtx_lock(&ccLock); 1540 1541 if(*cache == NULL) { 1542 createCache(cache, status); 1543 if(U_FAILURE(status)) { 1544 umtx_unlock(&ccLock); 1545 return 0; 1546 } 1547 } 1548 1549 res = uhash_igeti((*cache)->fTable, key); 1550 U_DEBUG_ASTRO_MSG(("%p: GET: [%d] == %d\n", (*cache)->fTable, key, res)); 1551 1552 umtx_unlock(&ccLock); 1553 return res; 1554 } 1555 1556 void CalendarCache::put(CalendarCache** cache, int32_t key, int32_t value, UErrorCode &status) { 1557 if(U_FAILURE(status)) { 1558 return; 1559 } 1560 umtx_lock(&ccLock); 1561 1562 if(*cache == NULL) { 1563 createCache(cache, status); 1564 if(U_FAILURE(status)) { 1565 umtx_unlock(&ccLock); 1566 return; 1567 } 1568 } 1569 1570 uhash_iputi((*cache)->fTable, key, value, &status); 1571 U_DEBUG_ASTRO_MSG(("%p: PUT: [%d] := %d\n", (*cache)->fTable, key, value)); 1572 1573 umtx_unlock(&ccLock); 1574 } 1575 1576 CalendarCache::CalendarCache(int32_t size, UErrorCode &status) { 1577 fTable = uhash_openSize(uhash_hashLong, uhash_compareLong, NULL, size, &status); 1578 U_DEBUG_ASTRO_MSG(("%p: Opening.\n", fTable)); 1579 } 1580 1581 CalendarCache::~CalendarCache() { 1582 if(fTable != NULL) { 1583 U_DEBUG_ASTRO_MSG(("%p: Closing.\n", fTable)); 1584 uhash_close(fTable); 1585 } 1586 } 1587 1588 U_NAMESPACE_END 1589 1590 #endif // !UCONFIG_NO_FORMATTING 1591
