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