jars/icu4j-4_2_1-src/src/com/ibm/icu/impl/CalendarAstronomer.java

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