sunposition.java

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

// **********************************************************************
// 
// <copyright>
// 
//  BBN Technologies
//  10 Moulton Street
//  Cambridge, MA 02138
//  [phone]
// 
//  Copyright (C) BBNT Solutions LLC. All rights reserved.
// 
// </copyright>
// **********************************************************************
// 
// $Source: /cvs/distapps/openmap/src/openmap/com/bbn/openmap/layer/daynight/SunPosition.java,v $
// $RCSfile: SunPosition.java,v $
// $Revision: 1.3.2.1 $
// $Date: 2004/10/14 18:27:03 $
// $Author: dietrick $
// 
// **********************************************************************

package com.bbn.openmap.layer.daynight;

import java.util.Calendar;
import java.util.Date;
import java.util.GregorianCalendar;

import com.bbn.openmap.MoreMath;
import com.bbn.openmap.LatLonPoint;

/**
 * SunPosition calculates the latitude/longitude on the Earth that is
 * closest to the Sun, the point on the Earth where the sun is
 * directly overhead.
 * <P>
 * 
 * All of these calculations are based on an epoch, or a starting
 * point where the Sun's position is known. From the reference book,
 * the epoch is 1990 January 0.0.
 * <P>
 * 
 * Also, all of the equations, and understanding of this whole
 * algorithm, was gained from a great book - Practical Astronomy with
 * Your Calculator, by Peter Duffett-Smith, Third Edition, Cambridge
 * University Press 1988.
 */
public class SunPosition {

    /** Epoch Julian Date. */
    public final static double EPOCH_JULIAN_DATE = 2447891.5;
    /** Epoch start time in seconds. */
    public final static double EPOCH_TIME_SECS = 631065600;

    /**
     * Constant denoting the number of radians an object would travel
     * if it orbited around the earth in a day.
     */
    public static double ORBIT_RADS_PER_DAY = MoreMath.TWO_PI / 365.242191;

    /**
     * Ecliptic Longitude of earth at 1990 January epoch. From
     * Duffett-Smith, chapter 46, table 6. (279.403303 degrees
     * converted to radians).
     */
    public static final double ECLIPTIC_LONGITUDE_EPOCH = 4.87650757893409;
    /**
     * Variable notation of ECLIPTIC_LONGITUDE_EPOCH from
     * Duffett-Smith.
     */
    public static final double epsilon_g = ECLIPTIC_LONGITUDE_EPOCH;

    /**
     * Ecliptic Longitude of of perigee. From Duffett-Smith, chapter
     * 46, table 6. (282.768422 degrees converted to radians).
     */
    public static final double ECLIPTIC_LONGITUDE_PERIGEE = 4.935239985213178;
    /**
     * Variable notation of ECLIPTIC_LONGITUDE_PERIGEE from
     * Duffett-Smith.
     */
    public static final double omega_bar_g = ECLIPTIC_LONGITUDE_PERIGEE;

    /**
     * Eccentricity of orbit, from Duffett-Smith, chapter 46, table 6.
     */
    public static final double ECCENTRICITY = 0.016713;

    /**
     * MEAN_OBLIQUITY_OF_EPOCH gives the mean obliquity of the
     * ecliptic, which is the angle between the planes of the equator
     * and the ecliptic. Using the algorithm described in
     * Duffett-Smith, chapter 27, this is calculated for the 1990
     * January epoch to be .4091155 radians (23.440592 degrees).
     */
    public static final double MEAN_OBLIQUITY_OF_EPOCH = .4091155;

    // These parameters are used in the Moon position calculations.

    /**
     * Moon parameter, from Duffett-Smith, chapter 65, table 10. In
     * radians.
     */
    public static final double MOON_EPOCH_MEAN_LONGITUDE = 318.351648 * Math.PI / 180.0;
    /**
     * The algorithm representation for the moon
     * MOON_EPOCH_MEAN_LONGITUDE, "l".
     */
    public static final double el0 = MOON_EPOCH_MEAN_LONGITUDE;

    /**
     * Moon parameter, from Duffett-Smith, chapter 65, table 10. In
     * radians.
     */
    public static final double PERIGEE_EPOCH_MEAN_LONGITUDE = 36.340410 * Math.PI / 180.0;
    /**
     * The algorithm representation for the moon
     * PERIGEE_EPOCH_MEAN_LONGITUDE.
     */
    public static final double P0 = PERIGEE_EPOCH_MEAN_LONGITUDE;

    /**
     * Moon parameter, from Duffett-Smith, chapter 65, table 10. In
     * radians.
     */
    public static final double NODE_EPOCH_MEAN_LONGITUDE = 318.510107 * Math.PI / 180.0;
    /**
     * The algorithm representation for the moon
     * NODE_EPOCH_MEAN_LONGITUDE.
     */
    public static final double N0 = NODE_EPOCH_MEAN_LONGITUDE;

    /**
     * Moon parameter, from Duffett-Smith, chapter 65, table 10. In
     * radians.
     */
    public static final double MOON_ORBIT_INCLINATION = 5.145396 * Math.PI / 180.0;
    /**
     * The algorithm representation for the moon
     * MOON_ORBIT_INCLINATION, "i".
     */
    public static final double eye = MOON_ORBIT_INCLINATION;

    /** Moon parameter, from Duffett-Smith, chapter 65, table 10. */
    public static final double MOON_ECCENTRICITY = .054900;
    /** Moon parameter, from Duffett-Smith, chapter 65, table 10. */
    public static final double MAJOR_AXIS_MOON_ORBIT = 384401; //km
    /**
     * Moon parameter, from Duffett-Smith, chapter 65, table 10. In
     * radians.
     */
    public static final double MOON_ANGULAR_SIZE = .5181 * Math.PI / 180.0;
    /**
     * Moon parameter, from Duffett-Smith, chapter 65, table 10. In
     * radians.
     */
    public static final double MOON_PARALLAX = .9507 * Math.PI / 180.0;

    /**
     * Use Kepllers's equation to find the eccentric anomaly. From
     * Duffett-Smith, chapter 47.
     * 
     * @param M the angle that the Sun has moved since it passed
     *        through perigee.
     */
    public static double eccentricAnomaly(double M) {
        double delta;
        double E = M;
        while (true) {
            delta = E - (ECCENTRICITY * Math.sin(E)) - M;

            if (Math.abs(delta) <= 1E-10)
                break;
            E -= (delta / (1.0 - (ECCENTRICITY * Math.cos(E))));
        }
        return E;
    }

    /**
     * Calculate the mean anomaly of sun, in radians. From
     * Duffett-Smith, chapter 47.
     * 
     * @param daysSinceEpoch number of days since 1990 January epoch.
     */
    protected static double sunMeanAnomaly(double daysSinceEpoch) {

        double N = ORBIT_RADS_PER_DAY * daysSinceEpoch;
        N = N % MoreMath.TWO_PI;
        if (N < 0)
            N += MoreMath.TWO_PI;

        double M0 = N + epsilon_g - omega_bar_g;
        if (M0 < 0)
            M0 += MoreMath.TWO_PI;
        return M0;
    }

    /**
     * Calculate the ecliptic longitude of sun, in radians. From
     * Duffett-Smith, chapter 47.
     * 
     * @param M0 sun's mean anomaly, calculated for the requested time
     *        relative to the 1990 epoch.
     */
    protected static double sunEclipticLongitude(double M0) {
        double E = eccentricAnomaly(M0);
        double v = 2 * Math.atan(Math.sqrt((1 + ECCENTRICITY)
                / (1 - ECCENTRICITY))
                * Math.tan(E / 2.0));
        double ret = v + omega_bar_g;
        ret = adjustWithin2PI(ret);
        return ret;
    }

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

        return Math.atan2(Math.sin(lambda) * cos_e - Math.tan(beta) * sin_e,
                Math.cos(lambda));
    }

    /**