1 // © 2016 and later: Unicode, Inc. and others. 2 // License & terms of use: http://www.unicode.org/copyright.html#License 3 /* 4 ******************************************************************************* 5 * Copyright (C) 1996-2011, International Business Machines Corporation and * 6 * others. All Rights Reserved. * 7 ******************************************************************************* 8 */ 9 10 package com.ibm.icu.impl; 11 12 import java.util.Date; 13 import java.util.TimeZone; 14 15 /** 16 * <code>CalendarAstronomer</code> is a class that can perform the calculations to 17 * determine the positions of the sun and moon, the time of sunrise and 18 * sunset, and other astronomy-related data. The calculations it performs 19 * are in some cases quite complicated, and this utility class saves you 20 * the trouble of worrying about them. 21 * <p> 22 * The measurement of time is a very important part of astronomy. Because 23 * astronomical bodies are constantly in motion, observations are only valid 24 * at a given moment in time. Accordingly, each <code>CalendarAstronomer</code> 25 * object has a <code>time</code> property that determines the date 26 * and time for which its calculations are performed. You can set and 27 * retrieve this property with {@link #setDate setDate}, {@link #getDate getDate} 28 * and related methods. 29 * <p> 30 * Almost all of the calculations performed by this class, or by any 31 * astronomer, are approximations to various degrees of accuracy. The 32 * calculations in this class are mostly modelled after those described 33 * in the book 34 * <a href="http://www.amazon.com/exec/obidos/ISBN=0521356997" target="_top"> 35 * Practical Astronomy With Your Calculator</a>, by Peter J. 36 * Duffett-Smith, Cambridge University Press, 1990. This is an excellent 37 * book, and if you want a greater understanding of how these calculations 38 * are performed it a very good, readable starting point. 39 * <p> 40 * <strong>WARNING:</strong> This class is very early in its development, and 41 * it is highly likely that its API will change to some degree in the future. 42 * At the moment, it basically does just enough to support {@link com.ibm.icu.util.IslamicCalendar} 43 * and {@link com.ibm.icu.util.ChineseCalendar}. 44 * 45 * @author Laura Werner 46 * @author Alan Liu 47 * @internal 48 */ 49 public class CalendarAstronomer { 50 51 //------------------------------------------------------------------------- 52 // Astronomical constants 53 //------------------------------------------------------------------------- 54 55 /** 56 * The number of standard hours in one sidereal day. 57 * Approximately 24.93. 58 * @internal 59 */ 60 public static final double SIDEREAL_DAY = 23.93446960027; 61 62 /** 63 * The number of sidereal hours in one mean solar day. 64 * Approximately 24.07. 65 * @internal 66 */ 67 public static final double SOLAR_DAY = 24.065709816; 68 69 /** 70 * The average number of solar days from one new moon to the next. This is the time 71 * it takes for the moon to return the same ecliptic longitude as the sun. 72 * It is longer than the sidereal month because the sun's longitude increases 73 * during the year due to the revolution of the earth around the sun. 74 * Approximately 29.53. 75 * 76 * @see #SIDEREAL_MONTH 77 * @internal 78 */ 79 public static final double SYNODIC_MONTH = 29.530588853; 80 81 /** 82 * The average number of days it takes 83 * for the moon to return to the same ecliptic longitude relative to the 84 * stellar background. This is referred to as the sidereal month. 85 * It is shorter than the synodic month due to 86 * the revolution of the earth around the sun. 87 * Approximately 27.32. 88 * 89 * @see #SYNODIC_MONTH 90 * @internal 91 */ 92 public static final double SIDEREAL_MONTH = 27.32166; 93 94 /** 95 * The average number number of days between successive vernal equinoxes. 96 * Due to the precession of the earth's 97 * axis, this is not precisely the same as the sidereal year. 98 * Approximately 365.24 99 * 100 * @see #SIDEREAL_YEAR 101 * @internal 102 */ 103 public static final double TROPICAL_YEAR = 365.242191; 104 105 /** 106 * The average number of days it takes 107 * for the sun to return to the same position against the fixed stellar 108 * background. This is the duration of one orbit of the earth about the sun 109 * as it would appear to an outside observer. 110 * Due to the precession of the earth's 111 * axis, this is not precisely the same as the tropical year. 112 * Approximately 365.25. 113 * 114 * @see #TROPICAL_YEAR 115 * @internal 116 */ 117 public static final double SIDEREAL_YEAR = 365.25636; 118 119 //------------------------------------------------------------------------- 120 // Time-related constants 121 //------------------------------------------------------------------------- 122 123 /** 124 * The number of milliseconds in one second. 125 * @internal 126 */ 127 public static final int SECOND_MS = 1000; 128 129 /** 130 * The number of milliseconds in one minute. 131 * @internal 132 */ 133 public static final int MINUTE_MS = 60*SECOND_MS; 134 135 /** 136 * The number of milliseconds in one hour. 137 * @internal 138 */ 139 public static final int HOUR_MS = 60*MINUTE_MS; 140 141 /** 142 * The number of milliseconds in one day. 143 * @internal 144 */ 145 public static final long DAY_MS = 24*HOUR_MS; 146 147 /** 148 * The start of the julian day numbering scheme used by astronomers, which 149 * is 1/1/4713 BC (Julian), 12:00 GMT. This is given as the number of milliseconds 150 * since 1/1/1970 AD (Gregorian), a negative number. 151 * Note that julian day numbers and 152 * the Julian calendar are <em>not</em> the same thing. Also note that 153 * julian days start at <em>noon</em>, not midnight. 154 * @internal 155 */ 156 public static final long JULIAN_EPOCH_MS = -210866760000000L; 157 158 // static { 159 // Calendar cal = new GregorianCalendar(TimeZone.getTimeZone("GMT")); 160 // cal.clear(); 161 // cal.set(cal.ERA, 0); 162 // cal.set(cal.YEAR, 4713); 163 // cal.set(cal.MONTH, cal.JANUARY); 164 // cal.set(cal.DATE, 1); 165 // cal.set(cal.HOUR_OF_DAY, 12); 166 // System.out.println("1.5 Jan 4713 BC = " + cal.getTime().getTime()); 167 168 // cal.clear(); 169 // cal.set(cal.YEAR, 2000); 170 // cal.set(cal.MONTH, cal.JANUARY); 171 // cal.set(cal.DATE, 1); 172 // cal.add(cal.DATE, -1); 173 // System.out.println("0.0 Jan 2000 = " + cal.getTime().getTime()); 174 // } 175 176 /** 177 * Milliseconds value for 0.0 January 2000 AD. 178 */ 179 static final long EPOCH_2000_MS = 946598400000L; 180 181 //------------------------------------------------------------------------- 182 // Assorted private data used for conversions 183 //------------------------------------------------------------------------- 184 185 // My own copies of these so compilers are more likely to optimize them away 186 static private final double PI = 3.14159265358979323846; 187 static private final double PI2 = PI * 2.0; 188 189 static private final double RAD_HOUR = 12 / PI; // radians -> hours 190 static private final double DEG_RAD = PI / 180; // degrees -> radians 191 static private final double RAD_DEG = 180 / PI; // radians -> degrees 192 193 //------------------------------------------------------------------------- 194 // Constructors 195 //------------------------------------------------------------------------- 196 197 /** 198 * Construct a new <code>CalendarAstronomer</code> object that is initialized to 199 * the current date and time. 200 * @internal 201 */ CalendarAstronomer()202 public CalendarAstronomer() { 203 this(System.currentTimeMillis()); 204 } 205 206 /** 207 * Construct a new <code>CalendarAstronomer</code> object that is initialized to 208 * the specified date and time. 209 * @internal 210 */ CalendarAstronomer(Date d)211 public CalendarAstronomer(Date d) { 212 this(d.getTime()); 213 } 214 215 /** 216 * Construct a new <code>CalendarAstronomer</code> object that is initialized to 217 * the specified time. The time is expressed as a number of milliseconds since 218 * January 1, 1970 AD (Gregorian). 219 * 220 * @see java.util.Date#getTime() 221 * @internal 222 */ CalendarAstronomer(long aTime)223 public CalendarAstronomer(long aTime) { 224 time = aTime; 225 } 226 227 /** 228 * Construct a new <code>CalendarAstronomer</code> object with the given 229 * latitude and longitude. The object's time is set to the current 230 * date and time. 231 * <p> 232 * @param longitude The desired longitude, in <em>degrees</em> east of 233 * the Greenwich meridian. 234 * 235 * @param latitude The desired latitude, in <em>degrees</em>. Positive 236 * values signify North, negative South. 237 * 238 * @see java.util.Date#getTime() 239 * @internal 240 */ CalendarAstronomer(double longitude, double latitude)241 public CalendarAstronomer(double longitude, double latitude) { 242 this(); 243 fLongitude = normPI(longitude * DEG_RAD); 244 fLatitude = normPI(latitude * DEG_RAD); 245 fGmtOffset = (long)(fLongitude * 24 * HOUR_MS / PI2); 246 } 247 248 249 //------------------------------------------------------------------------- 250 // Time and date getters and setters 251 //------------------------------------------------------------------------- 252 253 /** 254 * Set the current date and time of this <code>CalendarAstronomer</code> object. All 255 * astronomical calculations are performed based on this time setting. 256 * 257 * @param aTime the date and time, expressed as the number of milliseconds since 258 * 1/1/1970 0:00 GMT (Gregorian). 259 * 260 * @see #setDate 261 * @see #getTime 262 * @internal 263 */ setTime(long aTime)264 public void setTime(long aTime) { 265 time = aTime; 266 clearCache(); 267 } 268 269 /** 270 * Set the current date and time of this <code>CalendarAstronomer</code> object. All 271 * astronomical calculations are performed based on this time setting. 272 * 273 * @param date the time and date, expressed as a <code>Date</code> object. 274 * 275 * @see #setTime 276 * @see #getDate 277 * @internal 278 */ setDate(Date date)279 public void setDate(Date date) { 280 setTime(date.getTime()); 281 } 282 283 /** 284 * Set the current date and time of this <code>CalendarAstronomer</code> object. All 285 * astronomical calculations are performed based on this time setting. 286 * 287 * @param jdn the desired time, expressed as a "julian day number", 288 * which is the number of elapsed days since 289 * 1/1/4713 BC (Julian), 12:00 GMT. Note that julian day 290 * numbers start at <em>noon</em>. To get the jdn for 291 * the corresponding midnight, subtract 0.5. 292 * 293 * @see #getJulianDay 294 * @see #JULIAN_EPOCH_MS 295 * @internal 296 */ setJulianDay(double jdn)297 public void setJulianDay(double jdn) { 298 time = (long)(jdn * DAY_MS) + JULIAN_EPOCH_MS; 299 clearCache(); 300 julianDay = jdn; 301 } 302 303 /** 304 * Get the current time of this <code>CalendarAstronomer</code> object, 305 * represented as the number of milliseconds since 306 * 1/1/1970 AD 0:00 GMT (Gregorian). 307 * 308 * @see #setTime 309 * @see #getDate 310 * @internal 311 */ getTime()312 public long getTime() { 313 return time; 314 } 315 316 /** 317 * Get the current time of this <code>CalendarAstronomer</code> object, 318 * represented as a <code>Date</code> object. 319 * 320 * @see #setDate 321 * @see #getTime 322 * @internal 323 */ getDate()324 public Date getDate() { 325 return new Date(time); 326 } 327 328 /** 329 * Get the current time of this <code>CalendarAstronomer</code> object, 330 * expressed as a "julian day number", which is the number of elapsed 331 * days since 1/1/4713 BC (Julian), 12:00 GMT. 332 * 333 * @see #setJulianDay 334 * @see #JULIAN_EPOCH_MS 335 * @internal 336 */ getJulianDay()337 public double getJulianDay() { 338 if (julianDay == INVALID) { 339 julianDay = (double)(time - JULIAN_EPOCH_MS) / (double)DAY_MS; 340 } 341 return julianDay; 342 } 343 344 /** 345 * Return this object's time expressed in julian centuries: 346 * the number of centuries after 1/1/1900 AD, 12:00 GMT 347 * 348 * @see #getJulianDay 349 * @internal 350 */ getJulianCentury()351 public double getJulianCentury() { 352 if (julianCentury == INVALID) { 353 julianCentury = (getJulianDay() - 2415020.0) / 36525; 354 } 355 return julianCentury; 356 } 357 358 /** 359 * Returns the current Greenwich sidereal time, measured in hours 360 * @internal 361 */ getGreenwichSidereal()362 public double getGreenwichSidereal() { 363 if (siderealTime == INVALID) { 364 // See page 86 of "Practial Astronomy with your Calculator", 365 // by Peter Duffet-Smith, for details on the algorithm. 366 367 double UT = normalize((double)time/HOUR_MS, 24); 368 369 siderealTime = normalize(getSiderealOffset() + UT*1.002737909, 24); 370 } 371 return siderealTime; 372 } 373 getSiderealOffset()374 private double getSiderealOffset() { 375 if (siderealT0 == INVALID) { 376 double JD = Math.floor(getJulianDay() - 0.5) + 0.5; 377 double S = JD - 2451545.0; 378 double T = S / 36525.0; 379 siderealT0 = normalize(6.697374558 + 2400.051336*T + 0.000025862*T*T, 24); 380 } 381 return siderealT0; 382 } 383 384 /** 385 * Returns the current local sidereal time, measured in hours 386 * @internal 387 */ getLocalSidereal()388 public double getLocalSidereal() { 389 return normalize(getGreenwichSidereal() + (double)fGmtOffset/HOUR_MS, 24); 390 } 391 392 /** 393 * Converts local sidereal time to Universal Time. 394 * 395 * @param lst The Local Sidereal Time, in hours since sidereal midnight 396 * on this object's current date. 397 * 398 * @return The corresponding Universal Time, in milliseconds since 399 * 1 Jan 1970, GMT. 400 */ lstToUT(double lst)401 private long lstToUT(double lst) { 402 // Convert to local mean time 403 double lt = normalize((lst - getSiderealOffset()) * 0.9972695663, 24); 404 405 // Then find local midnight on this day 406 long base = DAY_MS * ((time + fGmtOffset)/DAY_MS) - fGmtOffset; 407 408 //out(" lt =" + lt + " hours"); 409 //out(" base=" + new Date(base)); 410 411 return base + (long)(lt * HOUR_MS); 412 } 413 414 415 //------------------------------------------------------------------------- 416 // Coordinate transformations, all based on the current time of this object 417 //------------------------------------------------------------------------- 418 419 /** 420 * Convert from ecliptic to equatorial coordinates. 421 * 422 * @param ecliptic A point in the sky in ecliptic coordinates. 423 * @return The corresponding point in equatorial coordinates. 424 * @internal 425 */ eclipticToEquatorial(Ecliptic ecliptic)426 public final Equatorial eclipticToEquatorial(Ecliptic ecliptic) 427 { 428 return eclipticToEquatorial(ecliptic.longitude, ecliptic.latitude); 429 } 430 431 /** 432 * Convert from ecliptic to equatorial coordinates. 433 * 434 * @param eclipLong The ecliptic longitude 435 * @param eclipLat The ecliptic latitude 436 * 437 * @return The corresponding point in equatorial coordinates. 438 * @internal 439 */ eclipticToEquatorial(double eclipLong, double eclipLat)440 public final Equatorial eclipticToEquatorial(double eclipLong, double eclipLat) 441 { 442 // See page 42 of "Practial Astronomy with your Calculator", 443 // by Peter Duffet-Smith, for details on the algorithm. 444 445 double obliq = eclipticObliquity(); 446 double sinE = Math.sin(obliq); 447 double cosE = Math.cos(obliq); 448 449 double sinL = Math.sin(eclipLong); 450 double cosL = Math.cos(eclipLong); 451 452 double sinB = Math.sin(eclipLat); 453 double cosB = Math.cos(eclipLat); 454 double tanB = Math.tan(eclipLat); 455 456 return new Equatorial(Math.atan2(sinL*cosE - tanB*sinE, cosL), 457 Math.asin(sinB*cosE + cosB*sinE*sinL) ); 458 } 459 460 /** 461 * Convert from ecliptic longitude to equatorial coordinates. 462 * 463 * @param eclipLong The ecliptic longitude 464 * 465 * @return The corresponding point in equatorial coordinates. 466 * @internal 467 */ eclipticToEquatorial(double eclipLong)468 public final Equatorial eclipticToEquatorial(double eclipLong) 469 { 470 return eclipticToEquatorial(eclipLong, 0); // TODO: optimize 471 } 472 473 /** 474 * @internal 475 */ eclipticToHorizon(double eclipLong)476 public Horizon eclipticToHorizon(double eclipLong) 477 { 478 Equatorial equatorial = eclipticToEquatorial(eclipLong); 479 480 double H = getLocalSidereal()*PI/12 - equatorial.ascension; // Hour-angle 481 482 double sinH = Math.sin(H); 483 double cosH = Math.cos(H); 484 double sinD = Math.sin(equatorial.declination); 485 double cosD = Math.cos(equatorial.declination); 486 double sinL = Math.sin(fLatitude); 487 double cosL = Math.cos(fLatitude); 488 489 double altitude = Math.asin(sinD*sinL + cosD*cosL*cosH); 490 double azimuth = Math.atan2(-cosD*cosL*sinH, sinD - sinL * Math.sin(altitude)); 491 492 return new Horizon(azimuth, altitude); 493 } 494 495 496 //------------------------------------------------------------------------- 497 // The Sun 498 //------------------------------------------------------------------------- 499 500 // 501 // Parameters of the Sun's orbit as of the epoch Jan 0.0 1990 502 // Angles are in radians (after multiplying by PI/180) 503 // 504 static final double JD_EPOCH = 2447891.5; // Julian day of epoch 505 506 static final double SUN_ETA_G = 279.403303 * PI/180; // Ecliptic longitude at epoch 507 static final double SUN_OMEGA_G = 282.768422 * PI/180; // Ecliptic longitude of perigee 508 static final double SUN_E = 0.016713; // Eccentricity of orbit 509 //double sunR0 = 1.495585e8; // Semi-major axis in KM 510 //double sunTheta0 = 0.533128 * PI/180; // Angular diameter at R0 511 512 // The following three methods, which compute the sun parameters 513 // given above for an arbitrary epoch (whatever time the object is 514 // set to), make only a small difference as compared to using the 515 // above constants. E.g., Sunset times might differ by ~12 516 // seconds. Furthermore, the eta-g computation is befuddled by 517 // Duffet-Smith's incorrect coefficients (p.86). I've corrected 518 // the first-order coefficient but the others may be off too - no 519 // way of knowing without consulting another source. 520 521 // /** 522 // * Return the sun's ecliptic longitude at perigee for the current time. 523 // * See Duffett-Smith, p. 86. 524 // * @return radians 525 // */ 526 // private double getSunOmegaG() { 527 // double T = getJulianCentury(); 528 // return (281.2208444 + (1.719175 + 0.000452778*T)*T) * DEG_RAD; 529 // } 530 531 // /** 532 // * Return the sun's ecliptic longitude for the current time. 533 // * See Duffett-Smith, p. 86. 534 // * @return radians 535 // */ 536 // private double getSunEtaG() { 537 // double T = getJulianCentury(); 538 // //return (279.6966778 + (36000.76892 + 0.0003025*T)*T) * DEG_RAD; 539 // // 540 // // The above line is from Duffett-Smith, and yields manifestly wrong 541 // // results. The below constant is derived empirically to match the 542 // // constant he gives for the 1990 EPOCH. 543 // // 544 // return (279.6966778 + (-0.3262541582718024 + 0.0003025*T)*T) * DEG_RAD; 545 // } 546 547 // /** 548 // * Return the sun's eccentricity of orbit for the current time. 549 // * See Duffett-Smith, p. 86. 550 // * @return double 551 // */ 552 // private double getSunE() { 553 // double T = getJulianCentury(); 554 // return 0.01675104 - (0.0000418 + 0.000000126*T)*T; 555 // } 556 557 /** 558 * The longitude of the sun at the time specified by this object. 559 * The longitude is measured in radians along the ecliptic 560 * from the "first point of Aries," the point at which the ecliptic 561 * crosses the earth's equatorial plane at the vernal equinox. 562 * <p> 563 * Currently, this method uses an approximation of the two-body Kepler's 564 * equation for the earth and the sun. It does not take into account the 565 * perturbations caused by the other planets, the moon, etc. 566 * @internal 567 */ getSunLongitude()568 public double getSunLongitude() 569 { 570 // See page 86 of "Practial Astronomy with your Calculator", 571 // by Peter Duffet-Smith, for details on the algorithm. 572 573 if (sunLongitude == INVALID) { 574 double[] result = getSunLongitude(getJulianDay()); 575 sunLongitude = result[0]; 576 meanAnomalySun = result[1]; 577 } 578 return sunLongitude; 579 } 580 581 /** 582 * TODO Make this public when the entire class is package-private. 583 */ getSunLongitude(double julian)584 /*public*/ double[] getSunLongitude(double julian) 585 { 586 // See page 86 of "Practial Astronomy with your Calculator", 587 // by Peter Duffet-Smith, for details on the algorithm. 588 589 double day = julian - JD_EPOCH; // Days since epoch 590 591 // Find the angular distance the sun in a fictitious 592 // circular orbit has travelled since the epoch. 593 double epochAngle = norm2PI(PI2/TROPICAL_YEAR*day); 594 595 // The epoch wasn't at the sun's perigee; find the angular distance 596 // since perigee, which is called the "mean anomaly" 597 double meanAnomaly = norm2PI(epochAngle + SUN_ETA_G - SUN_OMEGA_G); 598 599 // Now find the "true anomaly", e.g. the real solar longitude 600 // by solving Kepler's equation for an elliptical orbit 601 // NOTE: The 3rd ed. of the book lists omega_g and eta_g in different 602 // equations; omega_g is to be correct. 603 return new double[] { 604 norm2PI(trueAnomaly(meanAnomaly, SUN_E) + SUN_OMEGA_G), 605 meanAnomaly 606 }; 607 } 608 609 /** 610 * The position of the sun at this object's current date and time, 611 * in equatorial coordinates. 612 * @internal 613 */ getSunPosition()614 public Equatorial getSunPosition() { 615 return eclipticToEquatorial(getSunLongitude(), 0); 616 } 617 618 private static class SolarLongitude { 619 double value; SolarLongitude(double val)620 SolarLongitude(double val) { value = val; } 621 } 622 623 /** 624 * Constant representing the vernal equinox. 625 * For use with {@link #getSunTime(SolarLongitude, boolean) getSunTime}. 626 * Note: In this case, "vernal" refers to the northern hemisphere's seasons. 627 * @internal 628 */ 629 public static final SolarLongitude VERNAL_EQUINOX = new SolarLongitude(0); 630 631 /** 632 * Constant representing the summer solstice. 633 * For use with {@link #getSunTime(SolarLongitude, boolean) getSunTime}. 634 * Note: In this case, "summer" refers to the northern hemisphere's seasons. 635 * @internal 636 */ 637 public static final SolarLongitude SUMMER_SOLSTICE = new SolarLongitude(PI/2); 638 639 /** 640 * Constant representing the autumnal equinox. 641 * For use with {@link #getSunTime(SolarLongitude, boolean) getSunTime}. 642 * Note: In this case, "autumn" refers to the northern hemisphere's seasons. 643 * @internal 644 */ 645 public static final SolarLongitude AUTUMN_EQUINOX = new SolarLongitude(PI); 646 647 /** 648 * Constant representing the winter solstice. 649 * For use with {@link #getSunTime(SolarLongitude, boolean) getSunTime}. 650 * Note: In this case, "winter" refers to the northern hemisphere's seasons. 651 * @internal 652 */ 653 public static final SolarLongitude WINTER_SOLSTICE = new SolarLongitude((PI*3)/2); 654 655 /** 656 * Find the next time at which the sun's ecliptic longitude will have 657 * the desired value. 658 * @internal 659 */ getSunTime(double desired, boolean next)660 public long getSunTime(double desired, boolean next) 661 { 662 return timeOfAngle( new AngleFunc() { @Override 663 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 * @internal 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 * @internal 690 */ 691 public long getSunRiseSet(boolean rise) { 692 long t0 = time; 693 694 // Make a rough guess: 6am or 6pm local time on the current day 695 long noon = ((time + fGmtOffset)/DAY_MS)*DAY_MS - fGmtOffset + 12*HOUR_MS; 696 697 setTime(noon + (rise ? -6L : 6L) * HOUR_MS); 698 699 long t = riseOrSet(new CoordFunc() { 700 @Override 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 * @internal 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 * @internal 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 * @internal 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 * @internal 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 * @internal 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 * @internal 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 * @internal 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 * @internal 1166 */ 1167 public long getMoonTime(double desired, boolean next) 1168 { 1169 return timeOfAngle( new AngleFunc() { 1170 @Override 1171 public double eval() { return getMoonAge(); } }, 1172 desired, 1173 SYNODIC_MONTH, 1174 MINUTE_MS, 1175 next); 1176 } 1177 1178 /** 1179 * Find the next or previous time at which the moon will be in the 1180 * desired phase. 1181 * <p> 1182 * @param desired The desired phase of the moon. 1183 * @param next <tt>true</tt> if the next occurrance of the phase 1184 * is desired, <tt>false</tt> for the previous occurrance. 1185 * @internal 1186 */ 1187 public long getMoonTime(MoonAge desired, boolean next) { 1188 return getMoonTime(desired.value, next); 1189 } 1190 1191 /** 1192 * Returns the time (GMT) of sunrise or sunset on the local date to which 1193 * this calendar is currently set. 1194 * @internal 1195 */ 1196 public long getMoonRiseSet(boolean rise) 1197 { 1198 return riseOrSet(new CoordFunc() { 1199 @Override 1200 public Equatorial eval() { return getMoonPosition(); } 1201 }, 1202 rise, 1203 .533 * DEG_RAD, // Angular Diameter 1204 34 /60.0 * DEG_RAD, // Refraction correction 1205 MINUTE_MS); // Desired accuracy 1206 } 1207 1208 //------------------------------------------------------------------------- 1209 // Interpolation methods for finding the time at which a given event occurs 1210 //------------------------------------------------------------------------- 1211 1212 private interface AngleFunc { 1213 public double eval(); 1214 } 1215 1216 private long timeOfAngle(AngleFunc func, double desired, 1217 double periodDays, long epsilon, boolean next) 1218 { 1219 // Find the value of the function at the current time 1220 double lastAngle = func.eval(); 1221 1222 // Find out how far we are from the desired angle 1223 double deltaAngle = norm2PI(desired - lastAngle) ; 1224 1225 // Using the average period, estimate the next (or previous) time at 1226 // which the desired angle occurs. 1227 double deltaT = (deltaAngle + (next ? 0 : -PI2)) * (periodDays*DAY_MS) / PI2; 1228 1229 double lastDeltaT = deltaT; // Liu 1230 long startTime = time; // Liu 1231 1232 setTime(time + (long)deltaT); 1233 1234 // Now iterate until we get the error below epsilon. Throughout 1235 // this loop we use normPI to get values in the range -Pi to Pi, 1236 // since we're using them as correction factors rather than absolute angles. 1237 do { 1238 // Evaluate the function at the time we've estimated 1239 double angle = func.eval(); 1240 1241 // Find the # of milliseconds per radian at this point on the curve 1242 double factor = Math.abs(deltaT / normPI(angle-lastAngle)); 1243 1244 // Correct the time estimate based on how far off the angle is 1245 deltaT = normPI(desired - angle) * factor; 1246 1247 // HACK: 1248 // 1249 // If abs(deltaT) begins to diverge we need to quit this loop. 1250 // This only appears to happen when attempting to locate, for 1251 // example, a new moon on the day of the new moon. E.g.: 1252 // 1253 // This result is correct: 1254 // newMoon(7508(Mon Jul 23 00:00:00 CST 1990,false))= 1255 // Sun Jul 22 10:57:41 CST 1990 1256 // 1257 // But attempting to make the same call a day earlier causes deltaT 1258 // to diverge: 1259 // CalendarAstronomer.timeOfAngle() diverging: 1.348508727575625E9 -> 1260 // 1.3649828540224032E9 1261 // newMoon(7507(Sun Jul 22 00:00:00 CST 1990,false))= 1262 // Sun Jul 08 13:56:15 CST 1990 1263 // 1264 // As a temporary solution, we catch this specific condition and 1265 // adjust our start time by one eighth period days (either forward 1266 // or backward) and try again. 1267 // Liu 11/9/00 1268 if (Math.abs(deltaT) > Math.abs(lastDeltaT)) { 1269 long delta = (long) (periodDays * DAY_MS / 8); 1270 setTime(startTime + (next ? delta : -delta)); 1271 return timeOfAngle(func, desired, periodDays, epsilon, next); 1272 } 1273 1274 lastDeltaT = deltaT; 1275 lastAngle = angle; 1276 1277 setTime(time + (long)deltaT); 1278 } 1279 while (Math.abs(deltaT) > epsilon); 1280 1281 return time; 1282 } 1283 1284 private interface CoordFunc { 1285 public Equatorial eval(); 1286 } 1287 1288 private long riseOrSet(CoordFunc func, boolean rise, 1289 double diameter, double refraction, 1290 long epsilon) 1291 { 1292 Equatorial pos = null; 1293 double tanL = Math.tan(fLatitude); 1294 long deltaT = Long.MAX_VALUE; 1295 int count = 0; 1296 1297 // 1298 // Calculate the object's position at the current time, then use that 1299 // position to calculate the time of rising or setting. The position 1300 // will be different at that time, so iterate until the error is allowable. 1301 // 1302 do { 1303 // See "Practical Astronomy With Your Calculator, section 33. 1304 pos = func.eval(); 1305 double angle = Math.acos(-tanL * Math.tan(pos.declination)); 1306 double lst = ((rise ? PI2-angle : angle) + pos.ascension ) * 24 / PI2; 1307 1308 // Convert from LST to Universal Time. 1309 long newTime = lstToUT( lst ); 1310 1311 deltaT = newTime - time; 1312 setTime(newTime); 1313 } 1314 while (++ count < 5 && Math.abs(deltaT) > epsilon); 1315 1316 // Calculate the correction due to refraction and the object's angular diameter 1317 double cosD = Math.cos(pos.declination); 1318 double psi = Math.acos(Math.sin(fLatitude) / cosD); 1319 double x = diameter / 2 + refraction; 1320 double y = Math.asin(Math.sin(x) / Math.sin(psi)); 1321 long delta = (long)((240 * y * RAD_DEG / cosD)*SECOND_MS); 1322 1323 return time + (rise ? -delta : delta); 1324 } 1325 1326 //------------------------------------------------------------------------- 1327 // Other utility methods 1328 //------------------------------------------------------------------------- 1329 1330 /*** 1331 * Given 'value', add or subtract 'range' until 0 <= 'value' < range. 1332 * The modulus operator. 1333 */ 1334 private static final double normalize(double value, double range) { 1335 return value - range * Math.floor(value / range); 1336 } 1337 1338 /** 1339 * Normalize an angle so that it's in the range 0 - 2pi. 1340 * For positive angles this is just (angle % 2pi), but the Java 1341 * mod operator doesn't work that way for negative numbers.... 1342 */ 1343 private static final double norm2PI(double angle) { 1344 return normalize(angle, PI2); 1345 } 1346 1347 /** 1348 * Normalize an angle into the range -PI - PI 1349 */ 1350 private static final double normPI(double angle) { 1351 return normalize(angle + PI, PI2) - PI; 1352 } 1353 1354 /** 1355 * Find the "true anomaly" (longitude) of an object from 1356 * its mean anomaly and the eccentricity of its orbit. This uses 1357 * an iterative solution to Kepler's equation. 1358 * 1359 * @param meanAnomaly The object's longitude calculated as if it were in 1360 * a regular, circular orbit, measured in radians 1361 * from the point of perigee. 1362 * 1363 * @param eccentricity The eccentricity of the orbit 1364 * 1365 * @return The true anomaly (longitude) measured in radians 1366 */ 1367 private double trueAnomaly(double meanAnomaly, double eccentricity) 1368 { 1369 // First, solve Kepler's equation iteratively 1370 // Duffett-Smith, p.90 1371 double delta; 1372 double E = meanAnomaly; 1373 do { 1374 delta = E - eccentricity * Math.sin(E) - meanAnomaly; 1375 E = E - delta / (1 - eccentricity * Math.cos(E)); 1376 } 1377 while (Math.abs(delta) > 1e-5); // epsilon = 1e-5 rad 1378 1379 return 2.0 * Math.atan( Math.tan(E/2) * Math.sqrt( (1+eccentricity) 1380 /(1-eccentricity) ) ); 1381 } 1382 1383 /** 1384 * Return the obliquity of the ecliptic (the angle between the ecliptic 1385 * and the earth's equator) at the current time. This varies due to 1386 * the precession of the earth's axis. 1387 * 1388 * @return the obliquity of the ecliptic relative to the equator, 1389 * measured in radians. 1390 */ 1391 private double eclipticObliquity() { 1392 if (eclipObliquity == INVALID) { 1393 final double epoch = 2451545.0; // 2000 AD, January 1.5 1394 1395 double T = (getJulianDay() - epoch) / 36525; 1396 1397 eclipObliquity = 23.439292 1398 - 46.815/3600 * T 1399 - 0.0006/3600 * T*T 1400 + 0.00181/3600 * T*T*T; 1401 1402 eclipObliquity *= DEG_RAD; 1403 } 1404 return eclipObliquity; 1405 } 1406 1407 1408 //------------------------------------------------------------------------- 1409 // Private data 1410 //------------------------------------------------------------------------- 1411 1412 /** 1413 * Current time in milliseconds since 1/1/1970 AD 1414 * @see java.util.Date#getTime 1415 */ 1416 private long time; 1417 1418 /* These aren't used yet, but they'll be needed for sunset calculations 1419 * and equatorial to horizon coordinate conversions 1420 */ 1421 private double fLongitude = 0.0; 1422 private double fLatitude = 0.0; 1423 private long fGmtOffset = 0; 1424 1425 // 1426 // The following fields are used to cache calculated results for improved 1427 // performance. These values all depend on the current time setting 1428 // of this object, so the clearCache method is provided. 1429 // 1430 static final private double INVALID = Double.MIN_VALUE; 1431 1432 private transient double julianDay = INVALID; 1433 private transient double julianCentury = INVALID; 1434 private transient double sunLongitude = INVALID; 1435 private transient double meanAnomalySun = INVALID; 1436 private transient double moonLongitude = INVALID; 1437 private transient double moonEclipLong = INVALID; 1438 //private transient double meanAnomalyMoon = INVALID; 1439 private transient double eclipObliquity = INVALID; 1440 private transient double siderealT0 = INVALID; 1441 private transient double siderealTime = INVALID; 1442 1443 private transient Equatorial moonPosition = null; 1444 1445 private void clearCache() { 1446 julianDay = INVALID; 1447 julianCentury = INVALID; 1448 sunLongitude = INVALID; 1449 meanAnomalySun = INVALID; 1450 moonLongitude = INVALID; 1451 moonEclipLong = INVALID; 1452 //meanAnomalyMoon = INVALID; 1453 eclipObliquity = INVALID; 1454 siderealTime = INVALID; 1455 siderealT0 = INVALID; 1456 moonPosition = null; 1457 } 1458 1459 //private static void out(String s) { 1460 // System.out.println(s); 1461 //} 1462 1463 //private static String deg(double rad) { 1464 // return Double.toString(rad * RAD_DEG); 1465 //} 1466 1467 //private static String hours(long ms) { 1468 // return Double.toString((double)ms / HOUR_MS) + " hours"; 1469 //} 1470 1471 /** 1472 * @internal 1473 */ 1474 public String local(long localMillis) { 1475 return new Date(localMillis - TimeZone.getDefault().getRawOffset()).toString(); 1476 } 1477 1478 1479 /** 1480 * Represents the position of an object in the sky relative to the ecliptic, 1481 * the plane of the earth's orbit around the Sun. 1482 * This is a spherical coordinate system in which the latitude 1483 * specifies the position north or south of the plane of the ecliptic. 1484 * The longitude specifies the position along the ecliptic plane 1485 * relative to the "First Point of Aries", which is the Sun's position in the sky 1486 * at the Vernal Equinox. 1487 * <p> 1488 * Note that Ecliptic objects are immutable and cannot be modified 1489 * once they are constructed. This allows them to be passed and returned by 1490 * value without worrying about whether other code will modify them. 1491 * 1492 * @see CalendarAstronomer.Equatorial 1493 * @see CalendarAstronomer.Horizon 1494 * @internal 1495 */ 1496 public static final class Ecliptic { 1497 /** 1498 * Constructs an Ecliptic coordinate object. 1499 * <p> 1500 * @param lat The ecliptic latitude, measured in radians. 1501 * @param lon The ecliptic longitude, measured in radians. 1502 * @internal 1503 */ 1504 public Ecliptic(double lat, double lon) { 1505 latitude = lat; 1506 longitude = lon; 1507 } 1508 1509 /** 1510 * Return a string representation of this object 1511 * @internal 1512 */ 1513 @Override 1514 public String toString() { 1515 return Double.toString(longitude*RAD_DEG) + "," + (latitude*RAD_DEG); 1516 } 1517 1518 /** 1519 * The ecliptic latitude, in radians. This specifies an object's 1520 * position north or south of the plane of the ecliptic, 1521 * with positive angles representing north. 1522 * @internal 1523 */ 1524 public final double latitude; 1525 1526 /** 1527 * The ecliptic longitude, in radians. 1528 * This specifies an object's position along the ecliptic plane 1529 * relative to the "First Point of Aries", which is the Sun's position 1530 * in the sky at the Vernal Equinox, 1531 * with positive angles representing east. 1532 * <p> 1533 * A bit of trivia: the first point of Aries is currently in the 1534 * constellation Pisces, due to the precession of the earth's axis. 1535 * @internal 1536 */ 1537 public final double longitude; 1538 } 1539 1540 /** 1541 * Represents the position of an 1542 * object in the sky relative to the plane of the earth's equator. 1543 * The <i>Right Ascension</i> specifies the position east or west 1544 * along the equator, relative to the sun's position at the vernal 1545 * equinox. The <i>Declination</i> is the position north or south 1546 * of the equatorial plane. 1547 * <p> 1548 * Note that Equatorial objects are immutable and cannot be modified 1549 * once they are constructed. This allows them to be passed and returned by 1550 * value without worrying about whether other code will modify them. 1551 * 1552 * @see CalendarAstronomer.Ecliptic 1553 * @see CalendarAstronomer.Horizon 1554 * @internal 1555 */ 1556 public static final class Equatorial { 1557 /** 1558 * Constructs an Equatorial coordinate object. 1559 * <p> 1560 * @param asc The right ascension, measured in radians. 1561 * @param dec The declination, measured in radians. 1562 * @internal 1563 */ 1564 public Equatorial(double asc, double dec) { 1565 ascension = asc; 1566 declination = dec; 1567 } 1568 1569 /** 1570 * Return a string representation of this object, with the 1571 * angles measured in degrees. 1572 * @internal 1573 */ 1574 @Override 1575 public String toString() { 1576 return Double.toString(ascension*RAD_DEG) + "," + (declination*RAD_DEG); 1577 } 1578 1579 /** 1580 * Return a string representation of this object with the right ascension 1581 * measured in hours, minutes, and seconds. 1582 * @internal 1583 */ 1584 public String toHmsString() { 1585 return radToHms(ascension) + "," + radToDms(declination); 1586 } 1587 1588 /** 1589 * The right ascension, in radians. 1590 * This is the position east or west along the equator 1591 * relative to the sun's position at the vernal equinox, 1592 * with positive angles representing East. 1593 * @internal 1594 */ 1595 public final double ascension; 1596 1597 /** 1598 * The declination, in radians. 1599 * This is the position north or south of the equatorial plane, 1600 * with positive angles representing north. 1601 * @internal 1602 */ 1603 public final double declination; 1604 } 1605 1606 /** 1607 * Represents the position of an object in the sky relative to 1608 * the local horizon. 1609 * The <i>Altitude</i> represents the object's elevation above the horizon, 1610 * with objects below the horizon having a negative altitude. 1611 * The <i>Azimuth</i> is the geographic direction of the object from the 1612 * observer's position, with 0 representing north. The azimuth increases 1613 * clockwise from north. 1614 * <p> 1615 * Note that Horizon objects are immutable and cannot be modified 1616 * once they are constructed. This allows them to be passed and returned by 1617 * value without worrying about whether other code will modify them. 1618 * 1619 * @see CalendarAstronomer.Ecliptic 1620 * @see CalendarAstronomer.Equatorial 1621 * @internal 1622 */ 1623 public static final class Horizon { 1624 /** 1625 * Constructs a Horizon coordinate object. 1626 * <p> 1627 * @param alt The altitude, measured in radians above the horizon. 1628 * @param azim The azimuth, measured in radians clockwise from north. 1629 * @internal 1630 */ 1631 public Horizon(double alt, double azim) { 1632 altitude = alt; 1633 azimuth = azim; 1634 } 1635 1636 /** 1637 * Return a string representation of this object, with the 1638 * angles measured in degrees. 1639 * @internal 1640 */ 1641 @Override 1642 public String toString() { 1643 return Double.toString(altitude*RAD_DEG) + "," + (azimuth*RAD_DEG); 1644 } 1645 1646 /** 1647 * The object's altitude above the horizon, in radians. 1648 * @internal 1649 */ 1650 public final double altitude; 1651 1652 /** 1653 * The object's direction, in radians clockwise from north. 1654 * @internal 1655 */ 1656 public final double azimuth; 1657 } 1658 1659 static private String radToHms(double angle) { 1660 int hrs = (int) (angle*RAD_HOUR); 1661 int min = (int)((angle*RAD_HOUR - hrs) * 60); 1662 int sec = (int)((angle*RAD_HOUR - hrs - min/60.0) * 3600); 1663 1664 return Integer.toString(hrs) + "h" + min + "m" + sec + "s"; 1665 } 1666 1667 static private String radToDms(double angle) { 1668 int deg = (int) (angle*RAD_DEG); 1669 int min = (int)((angle*RAD_DEG - deg) * 60); 1670 int sec = (int)((angle*RAD_DEG - deg - min/60.0) * 3600); 1671 1672 return Integer.toString(deg) + "\u00b0" + min + "'" + sec + "\""; 1673 } 1674 } 1675