0001: /*
0002: *******************************************************************************
0003: * Copyright (C) 1996-2006, International Business Machines Corporation and *
0004: * others. All Rights Reserved. *
0005: *******************************************************************************
0006: */
0007: package com.ibm.icu.impl;
0008:
0009: import java.util.*;
0010:
0011: /**
0012: * <code>CalendarAstronomer</code> is a class that can perform the calculations to
0013: * determine the positions of the sun and moon, the time of sunrise and
0014: * sunset, and other astronomy-related data. The calculations it performs
0015: * are in some cases quite complicated, and this utility class saves you
0016: * the trouble of worrying about them.
0017: * <p>
0018: * The measurement of time is a very important part of astronomy. Because
0019: * astronomical bodies are constantly in motion, observations are only valid
0020: * at a given moment in time. Accordingly, each <code>CalendarAstronomer</code>
0021: * object has a <code>time</code> property that determines the date
0022: * and time for which its calculations are performed. You can set and
0023: * retrieve this property with {@link #setDate setDate}, {@link #getDate getDate}
0024: * and related methods.
0025: * <p>
0026: * Almost all of the calculations performed by this class, or by any
0027: * astronomer, are approximations to various degrees of accuracy. The
0028: * calculations in this class are mostly modelled after those described
0029: * in the book
0030: * <a href="http://www.amazon.com/exec/obidos/ISBN=0521356997" target="_top">
0031: * Practical Astronomy With Your Calculator</a>, by Peter J.
0032: * Duffett-Smith, Cambridge University Press, 1990. This is an excellent
0033: * book, and if you want a greater understanding of how these calculations
0034: * are performed it a very good, readable starting point.
0035: * <p>
0036: * <strong>WARNING:</strong> This class is very early in its development, and
0037: * it is highly likely that its API will change to some degree in the future.
0038: * At the moment, it basically does just enough to support {@link com.ibm.icu.util.IslamicCalendar}
0039: * and {@link com.ibm.icu.util.ChineseCalendar}.
0040: *
0041: * @author Laura Werner
0042: * @author Alan Liu
0043: * @internal
0044: */
0045: public class CalendarAstronomer {
0046:
0047: //-------------------------------------------------------------------------
0048: // Astronomical constants
0049: //-------------------------------------------------------------------------
0050:
0051: /**
0052: * The number of standard hours in one sidereal day.
0053: * Approximately 24.93.
0054: * @internal
0055: */
0056: public static final double SIDEREAL_DAY = 23.93446960027;
0057:
0058: /**
0059: * The number of sidereal hours in one mean solar day.
0060: * Approximately 24.07.
0061: * @internal
0062: */
0063: public static final double SOLAR_DAY = 24.065709816;
0064:
0065: /**
0066: * The average number of solar days from one new moon to the next. This is the time
0067: * it takes for the moon to return the same ecliptic longitude as the sun.
0068: * It is longer than the sidereal month because the sun's longitude increases
0069: * during the year due to the revolution of the earth around the sun.
0070: * Approximately 29.53.
0071: *
0072: * @see #SIDEREAL_MONTH
0073: * @internal
0074: */
0075: public static final double SYNODIC_MONTH = 29.530588853;
0076:
0077: /**
0078: * The average number of days it takes
0079: * for the moon to return to the same ecliptic longitude relative to the
0080: * stellar background. This is referred to as the sidereal month.
0081: * It is shorter than the synodic month due to
0082: * the revolution of the earth around the sun.
0083: * Approximately 27.32.
0084: *
0085: * @see #SYNODIC_MONTH
0086: * @internal
0087: */
0088: public static final double SIDEREAL_MONTH = 27.32166;
0089:
0090: /**
0091: * The average number number of days between successive vernal equinoxes.
0092: * Due to the precession of the earth's
0093: * axis, this is not precisely the same as the sidereal year.
0094: * Approximately 365.24
0095: *
0096: * @see #SIDEREAL_YEAR
0097: * @internal
0098: */
0099: public static final double TROPICAL_YEAR = 365.242191;
0100:
0101: /**
0102: * The average number of days it takes
0103: * for the sun to return to the same position against the fixed stellar
0104: * background. This is the duration of one orbit of the earth about the sun
0105: * as it would appear to an outside observer.
0106: * Due to the precession of the earth's
0107: * axis, this is not precisely the same as the tropical year.
0108: * Approximately 365.25.
0109: *
0110: * @see #TROPICAL_YEAR
0111: * @internal
0112: */
0113: public static final double SIDEREAL_YEAR = 365.25636;
0114:
0115: //-------------------------------------------------------------------------
0116: // Time-related constants
0117: //-------------------------------------------------------------------------
0118:
0119: /**
0120: * The number of milliseconds in one second.
0121: * @internal
0122: */
0123: public static final int SECOND_MS = 1000;
0124:
0125: /**
0126: * The number of milliseconds in one minute.
0127: * @internal
0128: */
0129: public static final int MINUTE_MS = 60 * SECOND_MS;
0130:
0131: /**
0132: * The number of milliseconds in one hour.
0133: * @internal
0134: */
0135: public static final int HOUR_MS = 60 * MINUTE_MS;
0136:
0137: /**
0138: * The number of milliseconds in one day.
0139: * @internal
0140: */
0141: public static final long DAY_MS = 24 * HOUR_MS;
0142:
0143: /**
0144: * The start of the julian day numbering scheme used by astronomers, which
0145: * is 1/1/4713 BC (Julian), 12:00 GMT. This is given as the number of milliseconds
0146: * since 1/1/1970 AD (Gregorian), a negative number.
0147: * Note that julian day numbers and
0148: * the Julian calendar are <em>not</em> the same thing. Also note that
0149: * julian days start at <em>noon</em>, not midnight.
0150: * @internal
0151: */
0152: public static final long JULIAN_EPOCH_MS = -210866760000000L;
0153:
0154: // static {
0155: // Calendar cal = new GregorianCalendar(TimeZone.getTimeZone("GMT"));
0156: // cal.clear();
0157: // cal.set(cal.ERA, 0);
0158: // cal.set(cal.YEAR, 4713);
0159: // cal.set(cal.MONTH, cal.JANUARY);
0160: // cal.set(cal.DATE, 1);
0161: // cal.set(cal.HOUR_OF_DAY, 12);
0162: // System.out.println("1.5 Jan 4713 BC = " + cal.getTime().getTime());
0163:
0164: // cal.clear();
0165: // cal.set(cal.YEAR, 2000);
0166: // cal.set(cal.MONTH, cal.JANUARY);
0167: // cal.set(cal.DATE, 1);
0168: // cal.add(cal.DATE, -1);
0169: // System.out.println("0.0 Jan 2000 = " + cal.getTime().getTime());
0170: // }
0171:
0172: /**
0173: * Milliseconds value for 0.0 January 2000 AD.
0174: */
0175: static final long EPOCH_2000_MS = 946598400000L;
0176:
0177: //-------------------------------------------------------------------------
0178: // Assorted private data used for conversions
0179: //-------------------------------------------------------------------------
0180:
0181: // My own copies of these so compilers are more likely to optimize them away
0182: static private final double PI = 3.14159265358979323846;
0183: static private final double PI2 = PI * 2.0;
0184:
0185: static private final double RAD_HOUR = 12 / PI; // radians -> hours
0186: static private final double DEG_RAD = PI / 180; // degrees -> radians
0187: static private final double RAD_DEG = 180 / PI; // radians -> degrees
0188:
0189: //-------------------------------------------------------------------------
0190: // Constructors
0191: //-------------------------------------------------------------------------
0192:
0193: /**
0194: * Construct a new <code>CalendarAstronomer</code> object that is initialized to
0195: * the current date and time.
0196: * @internal
0197: */
0198: public CalendarAstronomer() {
0199: this (System.currentTimeMillis());
0200: }
0201:
0202: /**
0203: * Construct a new <code>CalendarAstronomer</code> object that is initialized to
0204: * the specified date and time.
0205: * @internal
0206: */
0207: public CalendarAstronomer(Date d) {
0208: this (d.getTime());
0209: }
0210:
0211: /**
0212: * Construct a new <code>CalendarAstronomer</code> object that is initialized to
0213: * the specified time. The time is expressed as a number of milliseconds since
0214: * January 1, 1970 AD (Gregorian).
0215: *
0216: * @see java.util.Date#getTime()
0217: * @internal
0218: */
0219: public CalendarAstronomer(long aTime) {
0220: time = aTime;
0221: }
0222:
0223: /**
0224: * Construct a new <code>CalendarAstronomer</code> object with the given
0225: * latitude and longitude. The object's time is set to the current
0226: * date and time.
0227: * <p>
0228: * @param longitude The desired longitude, in <em>degrees</em> east of
0229: * the Greenwich meridian.
0230: *
0231: * @param latitude The desired latitude, in <em>degrees</em>. Positive
0232: * values signify North, negative South.
0233: *
0234: * @see java.util.Date#getTime()
0235: * @internal
0236: */
0237: public CalendarAstronomer(double longitude, double latitude) {
0238: this ();
0239: fLongitude = normPI(longitude * DEG_RAD);
0240: fLatitude = normPI(latitude * DEG_RAD);
0241: fGmtOffset = (long) (fLongitude * 24 * HOUR_MS / PI2);
0242: }
0243:
0244: //-------------------------------------------------------------------------
0245: // Time and date getters and setters
0246: //-------------------------------------------------------------------------
0247:
0248: /**
0249: * Set the current date and time of this <code>CalendarAstronomer</code> object. All
0250: * astronomical calculations are performed based on this time setting.
0251: *
0252: * @param aTime the date and time, expressed as the number of milliseconds since
0253: * 1/1/1970 0:00 GMT (Gregorian).
0254: *
0255: * @see #setDate
0256: * @see #getTime
0257: * @internal
0258: */
0259: public void setTime(long aTime) {
0260: time = aTime;
0261: clearCache();
0262: }
0263:
0264: /**
0265: * Set the current date and time of this <code>CalendarAstronomer</code> object. All
0266: * astronomical calculations are performed based on this time setting.
0267: *
0268: * @param date the time and date, expressed as a <code>Date</code> object.
0269: *
0270: * @see #setTime
0271: * @see #getDate
0272: * @internal
0273: */
0274: public void setDate(Date date) {
0275: setTime(date.getTime());
0276: }
0277:
0278: /**
0279: * Set the current date and time of this <code>CalendarAstronomer</code> object. All
0280: * astronomical calculations are performed based on this time setting.
0281: *
0282: * @param jdn the desired time, expressed as a "julian day number",
0283: * which is the number of elapsed days since
0284: * 1/1/4713 BC (Julian), 12:00 GMT. Note that julian day
0285: * numbers start at <em>noon</em>. To get the jdn for
0286: * the corresponding midnight, subtract 0.5.
0287: *
0288: * @see #getJulianDay
0289: * @see #JULIAN_EPOCH_MS
0290: * @internal
0291: */
0292: public void setJulianDay(double jdn) {
0293: time = (long) (jdn * DAY_MS) + JULIAN_EPOCH_MS;
0294: clearCache();
0295: julianDay = jdn;
0296: }
0297:
0298: /**
0299: * Get the current time of this <code>CalendarAstronomer</code> object,
0300: * represented as the number of milliseconds since
0301: * 1/1/1970 AD 0:00 GMT (Gregorian).
0302: *
0303: * @see #setTime
0304: * @see #getDate
0305: * @internal
0306: */
0307: public long getTime() {
0308: return time;
0309: }
0310:
0311: /**
0312: * Get the current time of this <code>CalendarAstronomer</code> object,
0313: * represented as a <code>Date</code> object.
0314: *
0315: * @see #setDate
0316: * @see #getTime
0317: * @internal
0318: */
0319: public Date getDate() {
0320: return new Date(time);
0321: }
0322:
0323: /**
0324: * Get the current time of this <code>CalendarAstronomer</code> object,
0325: * expressed as a "julian day number", which is the number of elapsed
0326: * days since 1/1/4713 BC (Julian), 12:00 GMT.
0327: *
0328: * @see #setJulianDay
0329: * @see #JULIAN_EPOCH_MS
0330: * @internal
0331: */
0332: public double getJulianDay() {
0333: if (julianDay == INVALID) {
0334: julianDay = (double) (time - JULIAN_EPOCH_MS)
0335: / (double) DAY_MS;
0336: }
0337: return julianDay;
0338: }
0339:
0340: /**
0341: * Return this object's time expressed in julian centuries:
0342: * the number of centuries after 1/1/1900 AD, 12:00 GMT
0343: *
0344: * @see #getJulianDay
0345: * @internal
0346: */
0347: public double getJulianCentury() {
0348: if (julianCentury == INVALID) {
0349: julianCentury = (getJulianDay() - 2415020.0) / 36525;
0350: }
0351: return julianCentury;
0352: }
0353:
0354: /**
0355: * Returns the current Greenwich sidereal time, measured in hours
0356: * @internal
0357: */
0358: public double getGreenwichSidereal() {
0359: if (siderealTime == INVALID) {
0360: // See page 86 of "Practial Astronomy with your Calculator",
0361: // by Peter Duffet-Smith, for details on the algorithm.
0362:
0363: double UT = normalize((double) time / HOUR_MS, 24);
0364:
0365: siderealTime = normalize(getSiderealOffset() + UT
0366: * 1.002737909, 24);
0367: }
0368: return siderealTime;
0369: }
0370:
0371: private double getSiderealOffset() {
0372: if (siderealT0 == INVALID) {
0373: double JD = Math.floor(getJulianDay() - 0.5) + 0.5;
0374: double S = JD - 2451545.0;
0375: double T = S / 36525.0;
0376: siderealT0 = normalize(6.697374558 + 2400.051336 * T
0377: + 0.000025862 * T * T, 24);
0378: }
0379: return siderealT0;
0380: }
0381:
0382: /**
0383: * Returns the current local sidereal time, measured in hours
0384: * @internal
0385: */
0386: public double getLocalSidereal() {
0387: return normalize(getGreenwichSidereal() + (double) fGmtOffset
0388: / HOUR_MS, 24);
0389: }
0390:
0391: /**
0392: * Converts local sidereal time to Universal Time.
0393: *
0394: * @param lst The Local Sidereal Time, in hours since sidereal midnight
0395: * on this object's current date.
0396: *
0397: * @return The corresponding Universal Time, in milliseconds since
0398: * 1 Jan 1970, GMT.
0399: */
0400: private long lstToUT(double lst) {
0401: // Convert to local mean time
0402: double lt = normalize(
0403: (lst - getSiderealOffset()) * 0.9972695663, 24);
0404:
0405: // Then find local midnight on this day
0406: long base = DAY_MS * ((time + fGmtOffset) / DAY_MS)
0407: - fGmtOffset;
0408:
0409: //out(" lt =" + lt + " hours");
0410: //out(" base=" + new Date(base));
0411:
0412: return base + (long) (lt * HOUR_MS);
0413: }
0414:
0415: //-------------------------------------------------------------------------
0416: // Coordinate transformations, all based on the current time of this object
0417: //-------------------------------------------------------------------------
0418:
0419: /**
0420: * Convert from ecliptic to equatorial coordinates.
0421: *
0422: * @param ecliptic A point in the sky in ecliptic coordinates.
0423: * @return The corresponding point in equatorial coordinates.
0424: * @internal
0425: */
0426: public final Equatorial eclipticToEquatorial(Ecliptic ecliptic) {
0427: return eclipticToEquatorial(ecliptic.longitude,
0428: ecliptic.latitude);
0429: }
0430:
0431: /**
0432: * Convert from ecliptic to equatorial coordinates.
0433: *
0434: * @param eclipLong The ecliptic longitude
0435: * @param eclipLat The ecliptic latitude
0436: *
0437: * @return The corresponding point in equatorial coordinates.
0438: * @internal
0439: */
0440: public final Equatorial eclipticToEquatorial(double eclipLong,
0441: double eclipLat) {
0442: // See page 42 of "Practial Astronomy with your Calculator",
0443: // by Peter Duffet-Smith, for details on the algorithm.
0444:
0445: double obliq = eclipticObliquity();
0446: double sinE = Math.sin(obliq);
0447: double cosE = Math.cos(obliq);
0448:
0449: double sinL = Math.sin(eclipLong);
0450: double cosL = Math.cos(eclipLong);
0451:
0452: double sinB = Math.sin(eclipLat);
0453: double cosB = Math.cos(eclipLat);
0454: double tanB = Math.tan(eclipLat);
0455:
0456: return new Equatorial(Math.atan2(sinL * cosE - tanB * sinE,
0457: cosL), Math.asin(sinB * cosE + cosB * sinE * sinL));
0458: }
0459:
0460: /**
0461: * Convert from ecliptic longitude to equatorial coordinates.
0462: *
0463: * @param eclipLong The ecliptic longitude
0464: *
0465: * @return The corresponding point in equatorial coordinates.
0466: * @internal
0467: */
0468: public final Equatorial eclipticToEquatorial(double eclipLong) {
0469: return eclipticToEquatorial(eclipLong, 0); // TODO: optimize
0470: }
0471:
0472: /**
0473: * @internal
0474: */
0475: public Horizon eclipticToHorizon(double eclipLong) {
0476: Equatorial equatorial = eclipticToEquatorial(eclipLong);
0477:
0478: double H = getLocalSidereal() * PI / 12 - equatorial.ascension; // Hour-angle
0479:
0480: double sinH = Math.sin(H);
0481: double cosH = Math.cos(H);
0482: double sinD = Math.sin(equatorial.declination);
0483: double cosD = Math.cos(equatorial.declination);
0484: double sinL = Math.sin(fLatitude);
0485: double cosL = Math.cos(fLatitude);
0486:
0487: double altitude = Math.asin(sinD * sinL + cosD * cosL * cosH);
0488: double azimuth = Math.atan2(-cosD * cosL * sinH, sinD - sinL
0489: * Math.sin(altitude));
0490:
0491: return new Horizon(azimuth, altitude);
0492: }
0493:
0494: //-------------------------------------------------------------------------
0495: // The Sun
0496: //-------------------------------------------------------------------------
0497:
0498: //
0499: // Parameters of the Sun's orbit as of the epoch Jan 0.0 1990
0500: // Angles are in radians (after multiplying by PI/180)
0501: //
0502: static final double JD_EPOCH = 2447891.5; // Julian day of epoch
0503:
0504: static final double SUN_ETA_G = 279.403303 * PI / 180; // Ecliptic longitude at epoch
0505: static final double SUN_OMEGA_G = 282.768422 * PI / 180; // Ecliptic longitude of perigee
0506: static final double SUN_E = 0.016713; // Eccentricity of orbit
0507:
0508: //double sunR0 = 1.495585e8; // Semi-major axis in KM
0509: //double sunTheta0 = 0.533128 * PI/180; // Angular diameter at R0
0510:
0511: // The following three methods, which compute the sun parameters
0512: // given above for an arbitrary epoch (whatever time the object is
0513: // set to), make only a small difference as compared to using the
0514: // above constants. E.g., Sunset times might differ by ~12
0515: // seconds. Furthermore, the eta-g computation is befuddled by
0516: // Duffet-Smith's incorrect coefficients (p.86). I've corrected
0517: // the first-order coefficient but the others may be off too - no
0518: // way of knowing without consulting another source.
0519:
0520: // /**
0521: // * Return the sun's ecliptic longitude at perigee for the current time.
0522: // * See Duffett-Smith, p. 86.
0523: // * @return radians
0524: // */
0525: // private double getSunOmegaG() {
0526: // double T = getJulianCentury();
0527: // return (281.2208444 + (1.719175 + 0.000452778*T)*T) * DEG_RAD;
0528: // }
0529:
0530: // /**
0531: // * Return the sun's ecliptic longitude for the current time.
0532: // * See Duffett-Smith, p. 86.
0533: // * @return radians
0534: // */
0535: // private double getSunEtaG() {
0536: // double T = getJulianCentury();
0537: // //return (279.6966778 + (36000.76892 + 0.0003025*T)*T) * DEG_RAD;
0538: // //
0539: // // The above line is from Duffett-Smith, and yields manifestly wrong
0540: // // results. The below constant is derived empirically to match the
0541: // // constant he gives for the 1990 EPOCH.
0542: // //
0543: // return (279.6966778 + (-0.3262541582718024 + 0.0003025*T)*T) * DEG_RAD;
0544: // }
0545:
0546: // /**
0547: // * Return the sun's eccentricity of orbit for the current time.
0548: // * See Duffett-Smith, p. 86.
0549: // * @return double
0550: // */
0551: // private double getSunE() {
0552: // double T = getJulianCentury();
0553: // return 0.01675104 - (0.0000418 + 0.000000126*T)*T;
0554: // }
0555:
0556: /**
0557: * The longitude of the sun at the time specified by this object.
0558: * The longitude is measured in radians along the ecliptic
0559: * from the "first point of Aries," the point at which the ecliptic
0560: * crosses the earth's equatorial plane at the vernal equinox.
0561: * <p>
0562: * Currently, this method uses an approximation of the two-body Kepler's
0563: * equation for the earth and the sun. It does not take into account the
0564: * perturbations caused by the other planets, the moon, etc.
0565: * @internal
0566: */
0567: public double getSunLongitude() {
0568: // See page 86 of "Practial Astronomy with your Calculator",
0569: // by Peter Duffet-Smith, for details on the algorithm.
0570:
0571: if (sunLongitude == INVALID) {
0572: double[] result = getSunLongitude(getJulianDay());
0573: sunLongitude = result[0];
0574: meanAnomalySun = result[1];
0575: }
0576: return sunLongitude;
0577: }
0578:
0579: /**
0580: * TODO Make this public when the entire class is package-private.
0581: */
0582: /*public*/double[] getSunLongitude(double julianDay) {
0583: // See page 86 of "Practial Astronomy with your Calculator",
0584: // by Peter Duffet-Smith, for details on the algorithm.
0585:
0586: double day = julianDay - JD_EPOCH; // Days since epoch
0587:
0588: // Find the angular distance the sun in a fictitious
0589: // circular orbit has travelled since the epoch.
0590: double epochAngle = norm2PI(PI2 / TROPICAL_YEAR * day);
0591:
0592: // The epoch wasn't at the sun's perigee; find the angular distance
0593: // since perigee, which is called the "mean anomaly"
0594: double meanAnomaly = norm2PI(epochAngle + SUN_ETA_G
0595: - SUN_OMEGA_G);
0596:
0597: // Now find the "true anomaly", e.g. the real solar longitude
0598: // by solving Kepler's equation for an elliptical orbit
0599: // NOTE: The 3rd ed. of the book lists omega_g and eta_g in different
0600: // equations; omega_g is to be correct.
0601: return new double[] {
0602: norm2PI(trueAnomaly(meanAnomaly, SUN_E) + SUN_OMEGA_G),
0603: meanAnomaly };
0604: }
0605:
0606: /**
0607: * The position of the sun at this object's current date and time,
0608: * in equatorial coordinates.
0609: * @internal
0610: */
0611: public Equatorial getSunPosition() {
0612: return eclipticToEquatorial(getSunLongitude(), 0);
0613: }
0614:
0615: private static class SolarLongitude {
0616: double value;
0617:
0618: SolarLongitude(double val) {
0619: value = val;
0620: }
0621: }
0622:
0623: /**
0624: * Constant representing the vernal equinox.
0625: * For use with {@link #getSunTime(SolarLongitude, boolean) getSunTime}.
0626: * Note: In this case, "vernal" refers to the northern hemisphere's seasons.
0627: * @internal
0628: */
0629: public static final SolarLongitude VERNAL_EQUINOX = new SolarLongitude(
0630: 0);
0631:
0632: /**
0633: * Constant representing the summer solstice.
0634: * For use with {@link #getSunTime(SolarLongitude, boolean) getSunTime}.
0635: * Note: In this case, "summer" refers to the northern hemisphere's seasons.
0636: * @internal
0637: */
0638: public static final SolarLongitude SUMMER_SOLSTICE = new SolarLongitude(
0639: PI / 2);
0640:
0641: /**
0642: * Constant representing the autumnal equinox.
0643: * For use with {@link #getSunTime(SolarLongitude, boolean) getSunTime}.
0644: * Note: In this case, "autumn" refers to the northern hemisphere's seasons.
0645: * @internal
0646: */
0647: public static final SolarLongitude AUTUMN_EQUINOX = new SolarLongitude(
0648: PI);
0649:
0650: /**
0651: * Constant representing the winter solstice.
0652: * For use with {@link #getSunTime(SolarLongitude, boolean) getSunTime}.
0653: * Note: In this case, "winter" refers to the northern hemisphere's seasons.
0654: * @internal
0655: */
0656: public static final SolarLongitude WINTER_SOLSTICE = new SolarLongitude(
0657: (PI * 3) / 2);
0658:
0659: /**
0660: * Find the next time at which the sun's ecliptic longitude will have
0661: * the desired value.
0662: * @internal
0663: */
0664: public long getSunTime(double desired, boolean next) {
0665: return timeOfAngle(new AngleFunc() {
0666: public double eval() {
0667: return getSunLongitude();
0668: }
0669: }, desired, TROPICAL_YEAR, MINUTE_MS, next);
0670: }
0671:
0672: /**
0673: * Find the next time at which the sun's ecliptic longitude will have
0674: * the desired value.
0675: * @internal
0676: */
0677: public long getSunTime(SolarLongitude desired, boolean next) {
0678: return getSunTime(desired.value, next);
0679: }
0680:
0681: /**
0682: * Returns the time (GMT) of sunrise or sunset on the local date to which
0683: * this calendar is currently set.
0684: *
0685: * NOTE: This method only works well if this object is set to a
0686: * time near local noon. Because of variations between the local
0687: * official time zone and the geographic longitude, the
0688: * computation can flop over into an adjacent day if this object
0689: * is set to a time near local midnight.
0690: *
0691: * @internal
0692: */
0693: public long getSunRiseSet(boolean rise) {
0694: long t0 = time;
0695:
0696: // Make a rough guess: 6am or 6pm local time on the current day
0697: long noon = ((time + fGmtOffset) / DAY_MS) * DAY_MS
0698: - fGmtOffset + 12 * HOUR_MS;
0699:
0700: setTime(noon + (long) ((rise ? -6 : 6) * HOUR_MS));
0701:
0702: long t = riseOrSet(new CoordFunc() {
0703: public Equatorial eval() {
0704: return getSunPosition();
0705: }
0706: }, rise, .533 * DEG_RAD, // Angular Diameter
0707: 34 / 60.0 * DEG_RAD, // Refraction correction
0708: MINUTE_MS / 12); // Desired accuracy
0709:
0710: setTime(t0);
0711: return t;
0712: }
0713:
0714: // Commented out - currently unused. ICU 2.6, Alan
0715: // //-------------------------------------------------------------------------
0716: // // Alternate Sun Rise/Set
0717: // // See Duffett-Smith p.93
0718: // //-------------------------------------------------------------------------
0719: //
0720: // // This yields worse results (as compared to USNO data) than getSunRiseSet().
0721: // /**
0722: // * TODO Make this public when the entire class is package-private.
0723: // */
0724: // /*public*/ long getSunRiseSet2(boolean rise) {
0725: // // 1. Calculate coordinates of the sun's center for midnight
0726: // double jd = Math.floor(getJulianDay() - 0.5) + 0.5;
0727: // double[] sl = getSunLongitude(jd);
0728: // double lambda1 = sl[0];
0729: // Equatorial pos1 = eclipticToEquatorial(lambda1, 0);
0730: //
0731: // // 2. Add ... to lambda to get position 24 hours later
0732: // double lambda2 = lambda1 + 0.985647*DEG_RAD;
0733: // Equatorial pos2 = eclipticToEquatorial(lambda2, 0);
0734: //
0735: // // 3. Calculate LSTs of rising and setting for these two positions
0736: // double tanL = Math.tan(fLatitude);
0737: // double H = Math.acos(-tanL * Math.tan(pos1.declination));
0738: // double lst1r = (PI2 + pos1.ascension - H) * 24 / PI2;
0739: // double lst1s = (pos1.ascension + H) * 24 / PI2;
0740: // H = Math.acos(-tanL * Math.tan(pos2.declination));
0741: // double lst2r = (PI2-H + pos2.ascension ) * 24 / PI2;
0742: // double lst2s = (H + pos2.ascension ) * 24 / PI2;
0743: // if (lst1r > 24) lst1r -= 24;
0744: // if (lst1s > 24) lst1s -= 24;
0745: // if (lst2r > 24) lst2r -= 24;
0746: // if (lst2s > 24) lst2s -= 24;
0747: //
0748: // // 4. Convert LSTs to GSTs. If GST1 > GST2, add 24 to GST2.
0749: // double gst1r = lstToGst(lst1r);
0750: // double gst1s = lstToGst(lst1s);
0751: // double gst2r = lstToGst(lst2r);
0752: // double gst2s = lstToGst(lst2s);
0753: // if (gst1r > gst2r) gst2r += 24;
0754: // if (gst1s > gst2s) gst2s += 24;
0755: //
0756: // // 5. Calculate GST at 0h UT of this date
0757: // double t00 = utToGst(0);
0758: //
0759: // // 6. Calculate GST at 0h on the observer's longitude
0760: // double offset = Math.round(fLongitude*12/PI); // p.95 step 6; he _rounds_ to nearest 15 deg.
0761: // double t00p = t00 - offset*1.002737909;
0762: // if (t00p < 0) t00p += 24; // do NOT normalize
0763: //
0764: // // 7. Adjust
0765: // if (gst1r < t00p) {
0766: // gst1r += 24;
0767: // gst2r += 24;
0768: // }
0769: // if (gst1s < t00p) {
0770: // gst1s += 24;
0771: // gst2s += 24;
0772: // }
0773: //
0774: // // 8.
0775: // double gstr = (24.07*gst1r-t00*(gst2r-gst1r))/(24.07+gst1r-gst2r);
0776: // double gsts = (24.07*gst1s-t00*(gst2s-gst1s))/(24.07+gst1s-gst2s);
0777: //
0778: // // 9. Correct for parallax, refraction, and sun's diameter
0779: // double dec = (pos1.declination + pos2.declination) / 2;
0780: // double psi = Math.acos(Math.sin(fLatitude) / Math.cos(dec));
0781: // double x = 0.830725 * DEG_RAD; // parallax+refraction+diameter
0782: // double y = Math.asin(Math.sin(x) / Math.sin(psi)) * RAD_DEG;
0783: // double delta_t = 240 * y / Math.cos(dec) / 3600; // hours
0784: //
0785: // // 10. Add correction to GSTs, subtract from GSTr
0786: // gstr -= delta_t;
0787: // gsts += delta_t;
0788: //
0789: // // 11. Convert GST to UT and then to local civil time
0790: // double ut = gstToUt(rise ? gstr : gsts);
0791: // //System.out.println((rise?"rise=":"set=") + ut + ", delta_t=" + delta_t);
0792: // long midnight = DAY_MS * (time / DAY_MS); // Find UT midnight on this day
0793: // return midnight + (long) (ut * 3600000);
0794: // }
0795:
0796: // Commented out - currently unused. ICU 2.6, Alan
0797: // /**
0798: // * Convert local sidereal time to Greenwich sidereal time.
0799: // * Section 15. Duffett-Smith p.21
0800: // * @param lst in hours (0..24)
0801: // * @return GST in hours (0..24)
0802: // */
0803: // double lstToGst(double lst) {
0804: // double delta = fLongitude * 24 / PI2;
0805: // return normalize(lst - delta, 24);
0806: // }
0807:
0808: // Commented out - currently unused. ICU 2.6, Alan
0809: // /**
0810: // * Convert UT to GST on this date.
0811: // * Section 12. Duffett-Smith p.17
0812: // * @param ut in hours
0813: // * @return GST in hours
0814: // */
0815: // double utToGst(double ut) {
0816: // return normalize(getT0() + ut*1.002737909, 24);
0817: // }
0818:
0819: // Commented out - currently unused. ICU 2.6, Alan
0820: // /**
0821: // * Convert GST to UT on this date.
0822: // * Section 13. Duffett-Smith p.18
0823: // * @param gst in hours
0824: // * @return UT in hours
0825: // */
0826: // double gstToUt(double gst) {
0827: // return normalize(gst - getT0(), 24) * 0.9972695663;
0828: // }
0829:
0830: // Commented out - currently unused. ICU 2.6, Alan
0831: // double getT0() {
0832: // // Common computation for UT <=> GST
0833: //
0834: // // Find JD for 0h UT
0835: // double jd = Math.floor(getJulianDay() - 0.5) + 0.5;
0836: //
0837: // double s = jd - 2451545.0;
0838: // double t = s / 36525.0;
0839: // double t0 = 6.697374558 + (2400.051336 + 0.000025862*t)*t;
0840: // return t0;
0841: // }
0842:
0843: // Commented out - currently unused. ICU 2.6, Alan
0844: // //-------------------------------------------------------------------------
0845: // // Alternate Sun Rise/Set
0846: // // See sci.astro FAQ
0847: // // http://www.faqs.org/faqs/astronomy/faq/part3/section-5.html
0848: // //-------------------------------------------------------------------------
0849: //
0850: // // Note: This method appears to produce inferior accuracy as
0851: // // compared to getSunRiseSet().
0852: //
0853: // /**
0854: // * TODO Make this public when the entire class is package-private.
0855: // */
0856: // /*public*/ long getSunRiseSet3(boolean rise) {
0857: //
0858: // // Compute day number for 0.0 Jan 2000 epoch
0859: // double d = (double)(time - EPOCH_2000_MS) / DAY_MS;
0860: //
0861: // // Now compute the Local Sidereal Time, LST:
0862: // //
0863: // double LST = 98.9818 + 0.985647352 * d + /*UT*15 + long*/
0864: // fLongitude*RAD_DEG;
0865: // //
0866: // // (east long. positive). Note that LST is here expressed in degrees,
0867: // // where 15 degrees corresponds to one hour. Since LST really is an angle,
0868: // // it's convenient to use one unit---degrees---throughout.
0869: //
0870: // // COMPUTING THE SUN'S POSITION
0871: // // ----------------------------
0872: // //
0873: // // To be able to compute the Sun's rise/set times, you need to be able to
0874: // // compute the Sun's position at any time. First compute the "day
0875: // // number" d as outlined above, for the desired moment. Next compute:
0876: // //
0877: // double oblecl = 23.4393 - 3.563E-7 * d;
0878: // //
0879: // double w = 282.9404 + 4.70935E-5 * d;
0880: // double M = 356.0470 + 0.9856002585 * d;
0881: // double e = 0.016709 - 1.151E-9 * d;
0882: // //
0883: // // This is the obliquity of the ecliptic, plus some of the elements of
0884: // // the Sun's apparent orbit (i.e., really the Earth's orbit): w =
0885: // // argument of perihelion, M = mean anomaly, e = eccentricity.
0886: // // Semi-major axis is here assumed to be exactly 1.0 (while not strictly
0887: // // true, this is still an accurate approximation). Next compute E, the
0888: // // eccentric anomaly:
0889: // //
0890: // double E = M + e*(180/PI) * Math.sin(M*DEG_RAD) * ( 1.0 + e*Math.cos(M*DEG_RAD) );
0891: // //
0892: // // where E and M are in degrees. This is it---no further iterations are
0893: // // needed because we know e has a sufficiently small value. Next compute
0894: // // the true anomaly, v, and the distance, r:
0895: // //
0896: // /* r * cos(v) = */ double A = Math.cos(E*DEG_RAD) - e;
0897: // /* r * sin(v) = */ double B = Math.sqrt(1 - e*e) * Math.sin(E*DEG_RAD);
0898: // //
0899: // // and
0900: // //
0901: // // r = sqrt( A*A + B*B )
0902: // double v = Math.atan2( B, A )*RAD_DEG;
0903: // //
0904: // // The Sun's true longitude, slon, can now be computed:
0905: // //
0906: // double slon = v + w;
0907: // //
0908: // // Since the Sun is always at the ecliptic (or at least very very close to
0909: // // it), we can use simplified formulae to convert slon (the Sun's ecliptic
0910: // // longitude) to sRA and sDec (the Sun's RA and Dec):
0911: // //
0912: // // sin(slon) * cos(oblecl)
0913: // // tan(sRA) = -------------------------
0914: // // cos(slon)
0915: // //
0916: // // sin(sDec) = sin(oblecl) * sin(slon)
0917: // //
0918: // // As was the case when computing az, the Azimuth, if possible use an
0919: // // atan2() function to compute sRA.
0920: //
0921: // double sRA = Math.atan2(Math.sin(slon*DEG_RAD) * Math.cos(oblecl*DEG_RAD), Math.cos(slon*DEG_RAD))*RAD_DEG;
0922: //
0923: // double sin_sDec = Math.sin(oblecl*DEG_RAD) * Math.sin(slon*DEG_RAD);
0924: // double sDec = Math.asin(sin_sDec)*RAD_DEG;
0925: //
0926: // // COMPUTING RISE AND SET TIMES
0927: // // ----------------------------
0928: // //
0929: // // To compute when an object rises or sets, you must compute when it
0930: // // passes the meridian and the HA of rise/set. Then the rise time is
0931: // // the meridian time minus HA for rise/set, and the set time is the
0932: // // meridian time plus the HA for rise/set.
0933: // //
0934: // // To find the meridian time, compute the Local Sidereal Time at 0h local
0935: // // time (or 0h UT if you prefer to work in UT) as outlined above---name
0936: // // that quantity LST0. The Meridian Time, MT, will now be:
0937: // //
0938: // // MT = RA - LST0
0939: // double MT = normalize(sRA - LST, 360);
0940: // //
0941: // // where "RA" is the object's Right Ascension (in degrees!). If negative,
0942: // // add 360 deg to MT. If the object is the Sun, leave the time as it is,
0943: // // but if it's stellar, multiply MT by 365.2422/366.2422, to convert from
0944: // // sidereal to solar time. Now, compute HA for rise/set, name that
0945: // // quantity HA0:
0946: // //
0947: // // sin(h0) - sin(lat) * sin(Dec)
0948: // // cos(HA0) = ---------------------------------
0949: // // cos(lat) * cos(Dec)
0950: // //
0951: // // where h0 is the altitude selected to represent rise/set. For a purely
0952: // // mathematical horizon, set h0 = 0 and simplify to:
0953: // //
0954: // // cos(HA0) = - tan(lat) * tan(Dec)
0955: // //
0956: // // If you want to account for refraction on the atmosphere, set h0 = -35/60
0957: // // degrees (-35 arc minutes), and if you want to compute the rise/set times
0958: // // for the Sun's upper limb, set h0 = -50/60 (-50 arc minutes).
0959: // //
0960: // double h0 = -50/60 * DEG_RAD;
0961: //
0962: // double HA0 = Math.acos(
0963: // (Math.sin(h0) - Math.sin(fLatitude) * sin_sDec) /
0964: // (Math.cos(fLatitude) * Math.cos(sDec*DEG_RAD)))*RAD_DEG;
0965: //
0966: // // When HA0 has been computed, leave it as it is for the Sun but multiply
0967: // // by 365.2422/366.2422 for stellar objects, to convert from sidereal to
0968: // // solar time. Finally compute:
0969: // //
0970: // // Rise time = MT - HA0
0971: // // Set time = MT + HA0
0972: // //
0973: // // convert the times from degrees to hours by dividing by 15.
0974: // //
0975: // // If you'd like to check that your calculations are accurate or just
0976: // // need a quick result, check the USNO's Sun or Moon Rise/Set Table,
0977: // // <URL:http://aa.usno.navy.mil/AA/data/docs/RS_OneYear.html>.
0978: //
0979: // double result = MT + (rise ? -HA0 : HA0); // in degrees
0980: //
0981: // // Find UT midnight on this day
0982: // long midnight = DAY_MS * (time / DAY_MS);
0983: //
0984: // return midnight + (long) (result * 3600000 / 15);
0985: // }
0986:
0987: //-------------------------------------------------------------------------
0988: // The Moon
0989: //-------------------------------------------------------------------------
0990:
0991: static final double moonL0 = 318.351648 * PI / 180; // Mean long. at epoch
0992: static final double moonP0 = 36.340410 * PI / 180; // Mean long. of perigee
0993: static final double moonN0 = 318.510107 * PI / 180; // Mean long. of node
0994: static final double moonI = 5.145366 * PI / 180; // Inclination of orbit
0995: static final double moonE = 0.054900; // Eccentricity of orbit
0996:
0997: // These aren't used right now
0998: static final double moonA = 3.84401e5; // semi-major axis (km)
0999: static final double moonT0 = 0.5181 * PI / 180; // Angular size at distance A
1000: static final double moonPi = 0.9507 * PI / 180; // Parallax at distance A
1001:
1002: /**
1003: * The position of the moon at the time set on this
1004: * object, in equatorial coordinates.
1005: * @internal
1006: */
1007: public Equatorial getMoonPosition() {
1008: //
1009: // See page 142 of "Practial Astronomy with your Calculator",
1010: // by Peter Duffet-Smith, for details on the algorithm.
1011: //
1012: if (moonPosition == null) {
1013: // Calculate the solar longitude. Has the side effect of
1014: // filling in "meanAnomalySun" as well.
1015: double sunLongitude = getSunLongitude();
1016:
1017: //
1018: // Find the # of days since the epoch of our orbital parameters.
1019: // TODO: Convert the time of day portion into ephemeris time
1020: //
1021: double day = getJulianDay() - JD_EPOCH; // Days since epoch
1022:
1023: // Calculate the mean longitude and anomaly of the moon, based on
1024: // a circular orbit. Similar to the corresponding solar calculation.
1025: double meanLongitude = norm2PI(13.1763966 * PI / 180 * day
1026: + moonL0);
1027: double meanAnomalyMoon = norm2PI(meanLongitude - 0.1114041
1028: * PI / 180 * day - moonP0);
1029:
1030: //
1031: // Calculate the following corrections:
1032: // Evection: the sun's gravity affects the moon's eccentricity
1033: // Annual Eqn: variation in the effect due to earth-sun distance
1034: // A3: correction factor (for ???)
1035: //
1036: double evection = 1.2739
1037: * PI
1038: / 180
1039: * Math.sin(2 * (meanLongitude - sunLongitude)
1040: - meanAnomalyMoon);
1041: double annual = 0.1858 * PI / 180
1042: * Math.sin(meanAnomalySun);
1043: double a3 = 0.3700 * PI / 180 * Math.sin(meanAnomalySun);
1044:
1045: meanAnomalyMoon += evection - annual - a3;
1046:
1047: //
1048: // More correction factors:
1049: // center equation of the center correction
1050: // a4 yet another error correction (???)
1051: //
1052: // TODO: Skip the equation of the center correction and solve Kepler's eqn?
1053: //
1054: double center = 6.2886 * PI / 180
1055: * Math.sin(meanAnomalyMoon);
1056: double a4 = 0.2140 * PI / 180
1057: * Math.sin(2 * meanAnomalyMoon);
1058:
1059: // Now find the moon's corrected longitude
1060: moonLongitude = meanLongitude + evection + center - annual
1061: + a4;
1062:
1063: //
1064: // And finally, find the variation, caused by the fact that the sun's
1065: // gravitational pull on the moon varies depending on which side of
1066: // the earth the moon is on
1067: //
1068: double variation = 0.6583 * PI / 180
1069: * Math.sin(2 * (moonLongitude - sunLongitude));
1070:
1071: moonLongitude += variation;
1072:
1073: //
1074: // What we've calculated so far is the moon's longitude in the plane
1075: // of its own orbit. Now map to the ecliptic to get the latitude
1076: // and longitude. First we need to find the longitude of the ascending
1077: // node, the position on the ecliptic where it is crossed by the moon's
1078: // orbit as it crosses from the southern to the northern hemisphere.
1079: //
1080: double nodeLongitude = norm2PI(moonN0 - 0.0529539 * PI
1081: / 180 * day);
1082:
1083: nodeLongitude -= 0.16 * PI / 180 * Math.sin(meanAnomalySun);
1084:
1085: double y = Math.sin(moonLongitude - nodeLongitude);
1086: double x = Math.cos(moonLongitude - nodeLongitude);
1087:
1088: moonEclipLong = Math.atan2(y * Math.cos(moonI), x)
1089: + nodeLongitude;
1090: double moonEclipLat = Math.asin(y * Math.sin(moonI));
1091:
1092: moonPosition = eclipticToEquatorial(moonEclipLong,
1093: moonEclipLat);
1094: }
1095: return moonPosition;
1096: }
1097:
1098: /**
1099: * The "age" of the moon at the time specified in this object.
1100: * This is really the angle between the
1101: * current ecliptic longitudes of the sun and the moon,
1102: * measured in radians.
1103: *
1104: * @see #getMoonPhase
1105: * @internal
1106: */
1107: public double getMoonAge() {
1108: // See page 147 of "Practial Astronomy with your Calculator",
1109: // by Peter Duffet-Smith, for details on the algorithm.
1110: //
1111: // Force the moon's position to be calculated. We're going to use
1112: // some the intermediate results cached during that calculation.
1113: //
1114: getMoonPosition();
1115:
1116: return norm2PI(moonEclipLong - sunLongitude);
1117: }
1118:
1119: /**
1120: * Calculate the phase of the moon at the time set in this object.
1121: * The returned phase is a <code>double</code> in the range
1122: * <code>0 <= phase < 1</code>, interpreted as follows:
1123: * <ul>
1124: * <li>0.00: New moon
1125: * <li>0.25: First quarter
1126: * <li>0.50: Full moon
1127: * <li>0.75: Last quarter
1128: * </ul>
1129: *
1130: * @see #getMoonAge
1131: * @internal
1132: */
1133: public double getMoonPhase() {
1134: // See page 147 of "Practial Astronomy with your Calculator",
1135: // by Peter Duffet-Smith, for details on the algorithm.
1136: return 0.5 * (1 - Math.cos(getMoonAge()));
1137: }
1138:
1139: private static class MoonAge {
1140: double value;
1141:
1142: MoonAge(double val) {
1143: value = val;
1144: }
1145: }
1146:
1147: /**
1148: * Constant representing a new moon.
1149: * For use with {@link #getMoonTime(MoonAge, boolean) getMoonTime}
1150: * @internal
1151: */
1152: public static final MoonAge NEW_MOON = new MoonAge(0);
1153:
1154: /**
1155: * Constant representing the moon's first quarter.
1156: * For use with {@link #getMoonTime(MoonAge, boolean) getMoonTime}
1157: * @internal
1158: */
1159: public static final MoonAge FIRST_QUARTER = new MoonAge(PI / 2);
1160:
1161: /**
1162: * Constant representing a full moon.
1163: * For use with {@link #getMoonTime(MoonAge, boolean) getMoonTime}
1164: * @internal
1165: */
1166: public static final MoonAge FULL_MOON = new MoonAge(PI);
1167:
1168: /**
1169: * Constant representing the moon's last quarter.
1170: * For use with {@link #getMoonTime(MoonAge, boolean) getMoonTime}
1171: * @internal
1172: */
1173: public static final MoonAge LAST_QUARTER = new MoonAge((PI * 3) / 2);
1174:
1175: /**
1176: * Find the next or previous time at which the Moon's ecliptic
1177: * longitude will have the desired value.
1178: * <p>
1179: * @param desired The desired longitude.
1180: * @param next <tt>true</tt> if the next occurrance of the phase
1181: * is desired, <tt>false</tt> for the previous occurrance.
1182: * @internal
1183: */
1184: public long getMoonTime(double desired, boolean next) {
1185: return timeOfAngle(new AngleFunc() {
1186: public double eval() {
1187: return getMoonAge();
1188: }
1189: }, desired, SYNODIC_MONTH, MINUTE_MS, next);
1190: }
1191:
1192: /**
1193: * Find the next or previous time at which the moon will be in the
1194: * desired phase.
1195: * <p>
1196: * @param desired The desired phase of the moon.
1197: * @param next <tt>true</tt> if the next occurrance of the phase
1198: * is desired, <tt>false</tt> for the previous occurrance.
1199: * @internal
1200: */
1201: public long getMoonTime(MoonAge desired, boolean next) {
1202: return getMoonTime(desired.value, next);
1203: }
1204:
1205: /**
1206: * Returns the time (GMT) of sunrise or sunset on the local date to which
1207: * this calendar is currently set.
1208: * @internal
1209: */
1210: public long getMoonRiseSet(boolean rise) {
1211: return riseOrSet(new CoordFunc() {
1212: public Equatorial eval() {
1213: return getMoonPosition();
1214: }
1215: }, rise, .533 * DEG_RAD, // Angular Diameter
1216: 34 / 60.0 * DEG_RAD, // Refraction correction
1217: MINUTE_MS); // Desired accuracy
1218: }
1219:
1220: //-------------------------------------------------------------------------
1221: // Interpolation methods for finding the time at which a given event occurs
1222: //-------------------------------------------------------------------------
1223:
1224: private interface AngleFunc {
1225: public double eval();
1226: }
1227:
1228: private long timeOfAngle(AngleFunc func, double desired,
1229: double periodDays, long epsilon, boolean next) {
1230: // Find the value of the function at the current time
1231: double lastAngle = func.eval();
1232:
1233: // Find out how far we are from the desired angle
1234: double deltaAngle = norm2PI(desired - lastAngle);
1235:
1236: // Using the average period, estimate the next (or previous) time at
1237: // which the desired angle occurs.
1238: double deltaT = (deltaAngle + (next ? 0 : -PI2))
1239: * (periodDays * DAY_MS) / PI2;
1240:
1241: double lastDeltaT = deltaT; // Liu
1242: long startTime = time; // Liu
1243:
1244: setTime(time + (long) deltaT);
1245:
1246: // Now iterate until we get the error below epsilon. Throughout
1247: // this loop we use normPI to get values in the range -Pi to Pi,
1248: // since we're using them as correction factors rather than absolute angles.
1249: do {
1250: // Evaluate the function at the time we've estimated
1251: double angle = func.eval();
1252:
1253: // Find the # of milliseconds per radian at this point on the curve
1254: double factor = Math
1255: .abs(deltaT / normPI(angle - lastAngle));
1256:
1257: // Correct the time estimate based on how far off the angle is
1258: deltaT = normPI(desired - angle) * factor;
1259:
1260: // HACK:
1261: //
1262: // If abs(deltaT) begins to diverge we need to quit this loop.
1263: // This only appears to happen when attempting to locate, for
1264: // example, a new moon on the day of the new moon. E.g.:
1265: //
1266: // This result is correct:
1267: // newMoon(7508(Mon Jul 23 00:00:00 CST 1990,false))=
1268: // Sun Jul 22 10:57:41 CST 1990
1269: //
1270: // But attempting to make the same call a day earlier causes deltaT
1271: // to diverge:
1272: // CalendarAstronomer.timeOfAngle() diverging: 1.348508727575625E9 ->
1273: // 1.3649828540224032E9
1274: // newMoon(7507(Sun Jul 22 00:00:00 CST 1990,false))=
1275: // Sun Jul 08 13:56:15 CST 1990
1276: //
1277: // As a temporary solution, we catch this specific condition and
1278: // adjust our start time by one eighth period days (either forward
1279: // or backward) and try again.
1280: // Liu 11/9/00
1281: if (Math.abs(deltaT) > Math.abs(lastDeltaT)) {
1282: long delta = (long) (periodDays * DAY_MS / 8);
1283: setTime(startTime + (next ? delta : -delta));
1284: return timeOfAngle(func, desired, periodDays, epsilon,
1285: next);
1286: }
1287:
1288: lastDeltaT = deltaT;
1289: lastAngle = angle;
1290:
1291: setTime(time + (long) deltaT);
1292: } while (Math.abs(deltaT) > epsilon);
1293:
1294: return time;
1295: }
1296:
1297: private interface CoordFunc {
1298: public Equatorial eval();
1299: }
1300:
1301: private long riseOrSet(CoordFunc func, boolean rise,
1302: double diameter, double refraction, long epsilon) {
1303: Equatorial pos = null;
1304: double tanL = Math.tan(fLatitude);
1305: long deltaT = Long.MAX_VALUE;
1306: int count = 0;
1307:
1308: //
1309: // Calculate the object's position at the current time, then use that
1310: // position to calculate the time of rising or setting. The position
1311: // will be different at that time, so iterate until the error is allowable.
1312: //
1313: do {
1314: // See "Practical Astronomy With Your Calculator, section 33.
1315: pos = func.eval();
1316: double angle = Math.acos(-tanL * Math.tan(pos.declination));
1317: double lst = ((rise ? PI2 - angle : angle) + pos.ascension)
1318: * 24 / PI2;
1319:
1320: // Convert from LST to Universal Time.
1321: long newTime = lstToUT(lst);
1322:
1323: deltaT = newTime - time;
1324: setTime(newTime);
1325: } while (++count < 5 && Math.abs(deltaT) > epsilon);
1326:
1327: // Calculate the correction due to refraction and the object's angular diameter
1328: double cosD = Math.cos(pos.declination);
1329: double psi = Math.acos(Math.sin(fLatitude) / cosD);
1330: double x = diameter / 2 + refraction;
1331: double y = Math.asin(Math.sin(x) / Math.sin(psi));
1332: long delta = (long) ((240 * y * RAD_DEG / cosD) * SECOND_MS);
1333:
1334: return time + (rise ? -delta : delta);
1335: }
1336:
1337: //-------------------------------------------------------------------------
1338: // Other utility methods
1339: //-------------------------------------------------------------------------
1340:
1341: /***
1342: * Given 'value', add or subtract 'range' until 0 <= 'value' < range.
1343: * The modulus operator.
1344: */
1345: private static final double normalize(double value, double range) {
1346: return value - range * Math.floor(value / range);
1347: }
1348:
1349: /**
1350: * Normalize an angle so that it's in the range 0 - 2pi.
1351: * For positive angles this is just (angle % 2pi), but the Java
1352: * mod operator doesn't work that way for negative numbers....
1353: */
1354: private static final double norm2PI(double angle) {
1355: return normalize(angle, PI2);
1356: }
1357:
1358: /**
1359: * Normalize an angle into the range -PI - PI
1360: */
1361: private static final double normPI(double angle) {
1362: return normalize(angle + PI, PI2) - PI;
1363: }
1364:
1365: /**
1366: * Find the "true anomaly" (longitude) of an object from
1367: * its mean anomaly and the eccentricity of its orbit. This uses
1368: * an iterative solution to Kepler's equation.
1369: *
1370: * @param meanAnomaly The object's longitude calculated as if it were in
1371: * a regular, circular orbit, measured in radians
1372: * from the point of perigee.
1373: *
1374: * @param eccentricity The eccentricity of the orbit
1375: *
1376: * @return The true anomaly (longitude) measured in radians
1377: */
1378: private double trueAnomaly(double meanAnomaly, double eccentricity) {
1379: // First, solve Kepler's equation iteratively
1380: // Duffett-Smith, p.90
1381: double delta;
1382: double E = meanAnomaly;
1383: do {
1384: delta = E - eccentricity * Math.sin(E) - meanAnomaly;
1385: E = E - delta / (1 - eccentricity * Math.cos(E));
1386: } while (Math.abs(delta) > 1e-5); // epsilon = 1e-5 rad
1387:
1388: return 2.0 * Math.atan(Math.tan(E / 2)
1389: * Math.sqrt((1 + eccentricity) / (1 - eccentricity)));
1390: }
1391:
1392: /**
1393: * Return the obliquity of the ecliptic (the angle between the ecliptic
1394: * and the earth's equator) at the current time. This varies due to
1395: * the precession of the earth's axis.
1396: *
1397: * @return the obliquity of the ecliptic relative to the equator,
1398: * measured in radians.
1399: */
1400: private double eclipticObliquity() {
1401: if (eclipObliquity == INVALID) {
1402: final double epoch = 2451545.0; // 2000 AD, January 1.5
1403:
1404: double T = (getJulianDay() - epoch) / 36525;
1405:
1406: eclipObliquity = 23.439292 - 46.815 / 3600 * T - 0.0006
1407: / 3600 * T * T + 0.00181 / 3600 * T * T * T;
1408:
1409: eclipObliquity *= DEG_RAD;
1410: }
1411: return eclipObliquity;
1412: }
1413:
1414: //-------------------------------------------------------------------------
1415: // Private data
1416: //-------------------------------------------------------------------------
1417:
1418: /**
1419: * Current time in milliseconds since 1/1/1970 AD
1420: * @see java.util.Date#getTime
1421: */
1422: private long time;
1423:
1424: /* These aren't used yet, but they'll be needed for sunset calculations
1425: * and equatorial to horizon coordinate conversions
1426: */
1427: private double fLongitude = 0.0;
1428: private double fLatitude = 0.0;
1429: private long fGmtOffset = 0;
1430:
1431: //
1432: // The following fields are used to cache calculated results for improved
1433: // performance. These values all depend on the current time setting
1434: // of this object, so the clearCache method is provided.
1435: //
1436: static final private double INVALID = Double.MIN_VALUE;
1437:
1438: private transient double julianDay = INVALID;
1439: private transient double julianCentury = INVALID;
1440: private transient double sunLongitude = INVALID;
1441: private transient double meanAnomalySun = INVALID;
1442: private transient double moonLongitude = INVALID;
1443: private transient double moonEclipLong = INVALID;
1444: private transient double meanAnomalyMoon = INVALID;
1445: private transient double eclipObliquity = INVALID;
1446: private transient double siderealT0 = INVALID;
1447: private transient double siderealTime = INVALID;
1448:
1449: private transient Equatorial moonPosition = null;
1450:
1451: private void clearCache() {
1452: julianDay = INVALID;
1453: julianCentury = INVALID;
1454: sunLongitude = INVALID;
1455: meanAnomalySun = INVALID;
1456: moonLongitude = INVALID;
1457: moonEclipLong = INVALID;
1458: meanAnomalyMoon = INVALID;
1459: eclipObliquity = INVALID;
1460: siderealTime = INVALID;
1461: siderealT0 = INVALID;
1462: moonPosition = null;
1463: }
1464:
1465: //private static void out(String s) {
1466: // System.out.println(s);
1467: //}
1468:
1469: //private static String deg(double rad) {
1470: // return Double.toString(rad * RAD_DEG);
1471: //}
1472:
1473: //private static String hours(long ms) {
1474: // return Double.toString((double)ms / HOUR_MS) + " hours";
1475: //}
1476:
1477: /**
1478: * @internal
1479: */
1480: public String local(long localMillis) {
1481: return new Date(localMillis
1482: - TimeZone.getDefault().getRawOffset()).toString();
1483: }
1484:
1485: /**
1486: * Represents the position of an object in the sky relative to the ecliptic,
1487: * the plane of the earth's orbit around the Sun.
1488: * This is a spherical coordinate system in which the latitude
1489: * specifies the position north or south of the plane of the ecliptic.
1490: * The longitude specifies the position along the ecliptic plane
1491: * relative to the "First Point of Aries", which is the Sun's position in the sky
1492: * at the Vernal Equinox.
1493: * <p>
1494: * Note that Ecliptic objects are immutable and cannot be modified
1495: * once they are constructed. This allows them to be passed and returned by
1496: * value without worrying about whether other code will modify them.
1497: *
1498: * @see CalendarAstronomer.Equatorial
1499: * @see CalendarAstronomer.Horizon
1500: * @internal
1501: */
1502: public static final class Ecliptic {
1503: /**
1504: * Constructs an Ecliptic coordinate object.
1505: * <p>
1506: * @param lat The ecliptic latitude, measured in radians.
1507: * @param lon The ecliptic longitude, measured in radians.
1508: * @internal
1509: */
1510: public Ecliptic(double lat, double lon) {
1511: latitude = lat;
1512: longitude = lon;
1513: }
1514:
1515: /**
1516: * Return a string representation of this object
1517: * @internal
1518: */
1519: public String toString() {
1520: return Double.toString(longitude * RAD_DEG) + ","
1521: + (latitude * RAD_DEG);
1522: }
1523:
1524: /**
1525: * The ecliptic latitude, in radians. This specifies an object's
1526: * position north or south of the plane of the ecliptic,
1527: * with positive angles representing north.
1528: * @internal
1529: */
1530: public final double latitude;
1531:
1532: /**
1533: * The ecliptic longitude, in radians.
1534: * This specifies an object's position along the ecliptic plane
1535: * relative to the "First Point of Aries", which is the Sun's position
1536: * in the sky at the Vernal Equinox,
1537: * with positive angles representing east.
1538: * <p>
1539: * A bit of trivia: the first point of Aries is currently in the
1540: * constellation Pisces, due to the precession of the earth's axis.
1541: * @internal
1542: */
1543: public final double longitude;
1544: }
1545:
1546: /**
1547: * Represents the position of an
1548: * object in the sky relative to the plane of the earth's equator.
1549: * The <i>Right Ascension</i> specifies the position east or west
1550: * along the equator, relative to the sun's position at the vernal
1551: * equinox. The <i>Declination</i> is the position north or south
1552: * of the equatorial plane.
1553: * <p>
1554: * Note that Equatorial objects are immutable and cannot be modified
1555: * once they are constructed. This allows them to be passed and returned by
1556: * value without worrying about whether other code will modify them.
1557: *
1558: * @see CalendarAstronomer.Ecliptic
1559: * @see CalendarAstronomer.Horizon
1560: * @internal
1561: */
1562: public static final class Equatorial {
1563: /**
1564: * Constructs an Equatorial coordinate object.
1565: * <p>
1566: * @param asc The right ascension, measured in radians.
1567: * @param dec The declination, measured in radians.
1568: * @internal
1569: */
1570: public Equatorial(double asc, double dec) {
1571: ascension = asc;
1572: declination = dec;
1573: }
1574:
1575: /**
1576: * Return a string representation of this object, with the
1577: * angles measured in degrees.
1578: * @internal
1579: */
1580: public String toString() {
1581: return Double.toString(ascension * RAD_DEG) + ","
1582: + (declination * RAD_DEG);
1583: }
1584:
1585: /**
1586: * Return a string representation of this object with the right ascension
1587: * measured in hours, minutes, and seconds.
1588: * @internal
1589: */
1590: public String toHmsString() {
1591: return radToHms(ascension) + "," + radToDms(declination);
1592: }
1593:
1594: /**
1595: * The right ascension, in radians.
1596: * This is the position east or west along the equator
1597: * relative to the sun's position at the vernal equinox,
1598: * with positive angles representing East.
1599: * @internal
1600: */
1601: public final double ascension;
1602:
1603: /**
1604: * The declination, in radians.
1605: * This is the position north or south of the equatorial plane,
1606: * with positive angles representing north.
1607: * @internal
1608: */
1609: public final double declination;
1610: }
1611:
1612: /**
1613: * Represents the position of an object in the sky relative to
1614: * the local horizon.
1615: * The <i>Altitude</i> represents the object's elevation above the horizon,
1616: * with objects below the horizon having a negative altitude.
1617: * The <i>Azimuth</i> is the geographic direction of the object from the
1618: * observer's position, with 0 representing north. The azimuth increases
1619: * clockwise from north.
1620: * <p>
1621: * Note that Horizon objects are immutable and cannot be modified
1622: * once they are constructed. This allows them to be passed and returned by
1623: * value without worrying about whether other code will modify them.
1624: *
1625: * @see CalendarAstronomer.Ecliptic
1626: * @see CalendarAstronomer.Equatorial
1627: * @internal
1628: */
1629: public static final class Horizon {
1630: /**
1631: * Constructs a Horizon coordinate object.
1632: * <p>
1633: * @param alt The altitude, measured in radians above the horizon.
1634: * @param azim The azimuth, measured in radians clockwise from north.
1635: * @internal
1636: */
1637: public Horizon(double alt, double azim) {
1638: altitude = alt;
1639: azimuth = azim;
1640: }
1641:
1642: /**
1643: * Return a string representation of this object, with the
1644: * angles measured in degrees.
1645: * @internal
1646: */
1647: public String toString() {
1648: return Double.toString(altitude * RAD_DEG) + ","
1649: + (azimuth * RAD_DEG);
1650: }
1651:
1652: /**
1653: * The object's altitude above the horizon, in radians.
1654: * @internal
1655: */
1656: public final double altitude;
1657:
1658: /**
1659: * The object's direction, in radians clockwise from north.
1660: * @internal
1661: */
1662: public final double azimuth;
1663: }
1664:
1665: static private String radToHms(double angle) {
1666: int hrs = (int) (angle * RAD_HOUR);
1667: int min = (int) ((angle * RAD_HOUR - hrs) * 60);
1668: int sec = (int) ((angle * RAD_HOUR - hrs - min / 60.0) * 3600);
1669:
1670: return Integer.toString(hrs) + "h" + min + "m" + sec + "s";
1671: }
1672:
1673: static private String radToDms(double angle) {
1674: int deg = (int) (angle * RAD_DEG);
1675: int min = (int) ((angle * RAD_DEG - deg) * 60);
1676: int sec = (int) ((angle * RAD_DEG - deg - min / 60.0) * 3600);
1677:
1678: return Integer.toString(deg) + "\u00b0" + min + "'" + sec
1679: + "\"";
1680: }
1681: }
|