jars/icu4j-4_4_2-src/main/classes/core/src/com/ibm/icu/impl/CalendarAstronomer.java

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