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