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 29static void debug_astro_loc(const char *f, int32_t l) 30{ 31 fprintf(stderr, "%s:%d: ", f, l); 32} 33 34static 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" 43static 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 62static inline UBool isINVALID(double d) { 63 return(uprv_isNaN(d)); 64} 65 66static UMTX ccLock = NULL; 67 68U_CDECL_BEGIN 69static UBool calendar_astro_cleanup(void) { 70 umtx_destroy(&ccLock); 71 return TRUE; 72} 73U_CDECL_END 74 75U_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 */ 104const 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 202const 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 */ 213inline 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 */ 222inline 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 */ 229inline 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 */ 243CalendarAstronomer::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 */ 254CalendarAstronomer::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 */ 273CalendarAstronomer::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 281CalendarAstronomer::~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 */ 301void 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 */ 322void 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 */ 338UDate 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 */ 352double 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 */ 367double 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 */ 379double 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 391double 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 */ 406double 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 */ 419double 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 */ 445CalendarAstronomer::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 */ 460CalendarAstronomer::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 */ 490CalendarAstronomer::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 */ 499CalendarAstronomer::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 */ 595static 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 */ 623double 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 */ 665CalendarAstronomer::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 */ 688double 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 */ 710double CalendarAstronomer::WINTER_SOLSTICE() { 711 return ((CalendarAstronomer::PI*3)/2); 712} 713 714CalendarAstronomer::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 */ 722class SunTimeAngleFunc : public CalendarAstronomer::AngleFunc { 723public: 724 virtual double eval(CalendarAstronomer& a) { return a.getSunLongitude(); } 725}; 726 727UDate 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 737CalendarAstronomer::CoordFunc::~CoordFunc() {} 738 739class RiseSetCoordFunc : public CalendarAstronomer::CoordFunc { 740public: 741 virtual void eval(CalendarAstronomer::Equatorial& result, CalendarAstronomer&a) { a.getSunPosition(result); } 742}; 743 744UDate 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 */ 1059const 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 */ 1149double 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 */ 1176double 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 */ 1188const 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 */ 1208const 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 1218class MoonTimeAngleFunc : public CalendarAstronomer::AngleFunc { 1219public: 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 */ 1237UDate 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 */ 1257UDate CalendarAstronomer::getMoonTime(const CalendarAstronomer::MoonAge& desired, UBool next) { 1258 return getMoonTime(desired.value, next); 1259} 1260 1261class MoonRiseSetCoordFunc : public CalendarAstronomer::CoordFunc { 1262public: 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 */ 1272UDate 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 1286UDate 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 1354UDate 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 */ 1403double 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//------------------------------------------------------------------------- 1423void 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 1467UnicodeString 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 1478UnicodeString 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 1490UnicodeString 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 1520void 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 1533int32_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 1556void 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 1576CalendarCache::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 1581CalendarCache::~CalendarCache() { 1582 if(fTable != NULL) { 1583 U_DEBUG_ASTRO_MSG(("%p: Closing.\n", fTable)); 1584 uhash_close(fTable); 1585 } 1586} 1587 1588U_NAMESPACE_END 1589 1590#endif // !UCONFIG_NO_FORMATTING 1591