sunposition.java

openmap java写的开源数字地图程序. 用applet实现,可以像google map 那样放大缩小地图.

     * Conversion from ecliptic to equatorial coordinates for
     * declination. From Duffett-Smith, chapter 27.
     * 
     * @param lambda ecliptic longitude
     * @param beta ecliptic latitude
     */
    protected static double eclipticToEquatorialDeclination(double lambda,
                                                            double beta) {
        double sin_e = Math.sin(MEAN_OBLIQUITY_OF_EPOCH);
        double cos_e = Math.cos(MEAN_OBLIQUITY_OF_EPOCH);

        return Math.asin(Math.sin(beta) * cos_e + Math.cos(beta) * sin_e
                * Math.sin(lambda));
    }

    /**
     * Given a date from a gregorian calendar, give back a julian
     * date. From Duffet-Smith, chapter 4.
     * 
     * @param cal Gregorian calendar for requested date.
     * @return julian date of request.
     */
    public static double calculateJulianDate(GregorianCalendar cal) {
        int year = cal.get(Calendar.YEAR);
        int month = cal.get(Calendar.MONTH);
        int day = cal.get(Calendar.DAY_OF_MONTH);

        // Algorithm expects that the January is = 1, which is
        // different from the Java representation
        month++;

        if ((month == 1) || (month == 2)) {
            year -= 1;
            month += 12;
        }

        int A = year / 100;
        int B = (int) (2 - A + (A / 4));
        int C = (int) (365.25 * (float) year);
        int D = (int) (30.6001 * (float) (month + 1));

        double julianDate = (double) (B + C + D + day) + 1720994.5;

        return julianDate;
    }

    /**
     * Calculate the greenwich sidereal time (GST). From
     * Duffett-Smith, chapter 12.
     * 
     * @param julianDate julian date of request
     * @param time calendar reflecting local time zone change to
     *        greenwich
     * @return GST relative to unix epoch.
     */
    public static double greenwichSiderealTime(double julianDate,
                                               GregorianCalendar time) {

        double T = (julianDate - 2451545.0) / 36525.0;
        double T0 = 6.697374558 + (T * (2400.051336 + (T + 2.5862E-5)));

        T0 = T0 % 24.0;
        if (T0 < 0) {
            T0 += 24;
        }

        double UT = time.get(Calendar.HOUR_OF_DAY)
                + (time.get(Calendar.MINUTE) + time.get(Calendar.SECOND) / 60.0)
                / 60.0;

        T0 += UT * 1.002737909;

        T0 = T0 % 24.0;
        if (T0 < 0) {
            T0 += 24.0;
        }

        return T0;
    }

    /**
     * Given the number of milliseconds since the unix epoch, compute
     * position on the earth (lat, lon) such that sun is directly
     * overhead. From Duffet-Smith, chapter 46-47.
     * 
     * @param mssue milliseconds since unix epoch
     * @return LatLonPoint of the point on the earth that is closest.
     */
    public static LatLonPoint sunPosition(long mssue) {

        // Set the date and clock, based on the millisecond count:
        GregorianCalendar cal = new GregorianCalendar();
        cal.setTime(new Date(mssue));

        double julianDate = calculateJulianDate(cal);

        // Need to correct time to GMT
        long gmtOffset = cal.get(Calendar.ZONE_OFFSET);
        // thanks to Erhard...
        long dstOffset = cal.get(Calendar.DST_OFFSET); // ins.
                                                       // 12.04.99
        cal.setTime(new Date(mssue - (gmtOffset + dstOffset))); // rep.
                                                                // 12.04.99

        double numDaysSinceEpoch = ((mssue / 1000) - EPOCH_TIME_SECS)
                / (24.0f * 3600.0f);

        // M0 - mean anomaly of the sun
        double M0 = sunMeanAnomaly(numDaysSinceEpoch);
        // lambda
        double sunLongitude = sunEclipticLongitude(M0);
        // alpha
        double sunAscension = eclipticToEquatorialAscension(sunLongitude, 0.0);
        // delta
        double sunDeclination = eclipticToEquatorialDeclination(sunLongitude,
                0.0);

        double tmpAscension = sunAscension - (MoreMath.TWO_PI / 24)
                * greenwichSiderealTime(julianDate, cal);

        return new LatLonPoint((float) (sunDeclination), (float) (tmpAscension), true);
    }

    /**
     * Given the number of milliseconds since the unix epoch, compute
     * position on the earth (lat, lon) such that moon is directly
     * overhead. From Duffet-Smith, chapter 65. Note: This is acting
     * like it works, but I don't have anything to test it against. No
     * promises.
     * 
     * @param mssue milliseconds since unix epoch
     * @return LatLonPoint of the point on the earth that is closest.
     */
    public static LatLonPoint moonPosition(long mssue) {

        // Set the date and clock, based on the millisecond count:
        GregorianCalendar cal = new GregorianCalendar();
        cal.setTime(new Date(mssue));

        double julianDate = calculateJulianDate(cal);

        // Need to correct time to GMT
        long gmtOffset = cal.get(Calendar.ZONE_OFFSET);
        cal.setTime(new Date(mssue - gmtOffset));

        // Step 1,2
        double numDaysSinceEpoch = ((mssue / 1000) - EPOCH_TIME_SECS)
                / (24.0f * 3600.0f);
        // Step 3
        // M0 - mean anomaly of the sun
        double M0 = sunMeanAnomaly(numDaysSinceEpoch);
        // lambda
        double sunLongitude = sunEclipticLongitude(M0);
        // Step 4
        double el = (13.1763966 * numDaysSinceEpoch * Math.PI / 180) + el0;
        el = adjustWithin2PI(el);
        // Step 5
        double Mm = el - (.1114041 * numDaysSinceEpoch * Math.PI / 180) - P0;
        Mm = adjustWithin2PI(Mm);
        // Step 6
        double N = N0 - (.0529539 * numDaysSinceEpoch * Math.PI / 180);
        N = adjustWithin2PI(N);
        // Step 7
        double C = el - sunLongitude;
        double Ev = 1.2739 * Math.sin(2 * C - Mm);
        // Step 8
        double Ae = .1858 * Math.sin(M0);
        double A3 = .37 * Math.sin(M0);
        // Step 9
        double Mmp = Mm + Ev - Ae - A3;
        // Step 10
        double Ec = 6.2886 * Math.sin(Mmp);
        //燬tep 11
        double A4 = 0.214 * Math.sin(2 * Mmp);
        //燬tep 12
        double elp = el + Ev + Ec - Ae + A4;
        //燬tep 13
        double V = .6583 * Math.sin(2 * (elp - sunLongitude));
        //燬tep 14
        double elpp = elp + V;
        //燬tep 15
        double Np = N - (.16 * Math.sin(M0));
        //燬tep 16
        double y = Math.sin(elpp - Np) * Math.cos(eye);
        //燬tep 17
        double x = Math.cos(elpp - Np);
        //燬tep 18
        double amb = Math.atan2(y, x);
        //燬tep 19
        double lambda_m = amb + Np;
        //燬tep 20
        double beta_m = Math.asin(Math.sin(elpp - Np) * Math.sin(eye));
        //燬tep 21
        // alpha
        double moonAscension = eclipticToEquatorialAscension(lambda_m, beta_m);
        // delta
        double moonDeclination = eclipticToEquatorialDeclination(lambda_m,
                beta_m);

        double tmpAscension = moonAscension - (MoreMath.TWO_PI / 24)
                * greenwichSiderealTime(julianDate, cal);

        return new LatLonPoint((float) (moonDeclination), (float) (tmpAscension), true);
    }

    /**
     * Little function that resets the input to be within 0 - 2*PI, by
     * adding or subtracting 2PI as needed.
     * 
     * @param num The number to be modified, if needed.
     */
    protected static double adjustWithin2PI(double num) {
        if (num < 0) {
            do
                num += MoreMath.TWO_PI;
            while (num < 0);
        } else if (num > MoreMath.TWO_PI) {
            do
                num -= MoreMath.TWO_PI;
            while (num > MoreMath.TWO_PI);
        }
        return num;
    }

    public static void main(String[] arg) {
        System.out.println("Sun is over "
                + SunPosition.sunPosition(System.currentTimeMillis()));
        System.out.println("Moon is over "
                + SunPosition.moonPosition(System.currentTimeMillis()));
    }
}