/*
 *    GeoTools - The Open Source Java GIS Toolkit
 *    http://geotools.org
 * 
 *    (C) 2001-2008, Open Source Geospatial Foundation (OSGeo)
 *   
 *    This library is free software; you can redistribute it and/or
 *    modify it under the terms of the GNU Lesser General Public
 *    License as published by the Free Software Foundation;
 *    version 2.1 of the License.
 *
 *    This library is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *    Lesser General Public License for more details.
 *
 *    NOTE: permission has been given to the JScience project (http://www.jscience.org)
 *          to distribute this file under BSD-like license.
 */
package edu.wisc.ssec.mcidasv.util;

// J2SE dependencies
import java.awt.geom.Point2D;
import java.text.DateFormat;
import java.text.ParseException;
import java.util.Date;
import java.util.TimeZone;


/**
 * Calcule la position du soleil relativement à la position de l'observateur.
 * Cette classe reçoit en entrés les coordonnées spatio-temporelles de
 * l'observateur, soit:
 *
 * <TABLE border='0'><TR><TD valign="top">
 * &nbsp;<BR>
 * <UL>
 *   <LI>La longitude (en degrées) de l'observateur;</LI>
 *   <LI>La latitude (en degrées) de l'observateur;</LI>
 *   <LI>La date et heure en heure universelle (GMT).</LI>
 * </UL>
 *
 * La position du soleil calculée en sortie comprend les valeurs suivantes:
 *
 * <UL>
 *   <LI>L'azimuth du soleil (en degrés dans le sens des aiguilles d'une montre depuis le nord);</LI>
 *   <LI>L'élévation du soleil (en degrés par rapport a l'horizon).</LI>
 * </UL>
 * </TD>
 *
 * <TD><img src="doc-files/CelestialSphere.png"></TD>
 * </TR></TABLE>
 *
 * Les algorithmes utilisés dans cette classe sont des adaptations des algorithmes
 * en javascript écrit par le "National Oceanic and Atmospheric Administration,
 * Surface Radiation Research Branch". L'application original est le
 *
 * <a href="http://www.srrb.noaa.gov/highlights/sunrise/azel.html">Solar Position Calculator</a>.
 *
 * <p>
 * The approximations used in these programs are very good for years between
 * 1800 and 2100. Results should still be sufficiently accurate for the range
 * from -1000 to 3000. Outside of this range, results will be given, but the
 * potential for error is higher.
 *
 * @since 2.1
 *
 *
 * @source $URL$
 * @version $Id$
 * @author Remi Eve
 * @author Martin Desruisseaux (IRD)
 */
public class SunRelativePosition {
    /**
     * Number of milliseconds in a day.
     */
    private static final int DAY_MILLIS = 24*60*60*1000;

    /**
     * Valeur affectée lorsque un resultat n'est pas calculable du
     * fait de la nuit. Cette valeur concerne les valeurs de sorties
     * {@link #elevation} et {@link #azimuth}.
     */
    private static final double DARK = Double.NaN;

    /**
     * {@linkplain #getElevation Elevation angle} of astronomical twilight, in degrees.
     * Astronomical twilight is the time of morning or evening when the sun is 18° below
     * the horizon (solar elevation angle of -18°).
     */
    public static final double ASTRONOMICAL_TWILIGHT = -18;

    /**
     * {@linkplain #getElevation Elevation angle} of nautical twilight, in degrees.
     * Nautical twilight is the time of morning or evening when the sun is 12° below
     * the horizon (solar elevation angle of -12°).
     */
    public static final double NAUTICAL_TWILIGHT = -12;

    /**
     * {@linkplain #getElevation Elevation angle} of civil twilight, in degrees. Civil
     * twilight is the time of morning or evening when the sun is 6° below the horizon
     * (solar elevation angle of -6°).
     */
    public static final double CIVIL_TWILIGHT = -6;

    /**
     * Sun's {@linkplain #getElevation elevation angle} at twilight, in degrees.
     * Common values are defined for the
     * {@linkplain #ASTRONOMICAL_TWILIGHT astronomical twilight} (-18°),
     * {@linkplain #NAUTICAL_TWILIGHT nautical twilight} (-12°) and
     * {@linkplain #CIVIL_TWILIGHT civil twilight} (-6°).
     * If no twilight are defined, then this value is {@linkplain Double#NaN NaN}.
     * The {@linkplain #getElevation elevation} and {@linkplain #getAzimuth azimuth} are
     * set to {@linkplain Double#NaN NaN} when the sun elevation is below the twilight
     * value (i.e. during night). The default value is {@link #CIVIL_TWILIGHT}.
     */
    private double twilight = CIVIL_TWILIGHT;

    /**
     * Heure à laquelle le soleil est au plus haut dans la journée en millisecondes
     * écoulées depuis le 1er janvier 1970.
     */
    private long noonTime;
    
    /**
     * Azimuth du soleil, en degrés dans le sens des
     * aiguilles d'une montre depuis le nord.
     */
    private double azimuth;

    /**
     * Elévation du soleil, en degrés par rapport a l'horizon.
     */
    private double elevation;

    /**
     * Geographic coordinate where current elevation and azimuth were computed.
     * Value are in degrees of longitude or latitude.
     */
    private double latitude, longitude;

    /**
     * Date and time when the current elevation and azimuth were computed.
     * Value is in milliseconds ellapsed since midnight UTC, January 1st, 1970.
     */
    private long time = System.currentTimeMillis();

    /**
     * {@code true} is the elevation and azimuth are computed, or {@code false}
     * if they need to be computed.  This flag is set to {@code false} when the date
     * and/or the coordinate change.
     */
    private boolean updated;

    /**
     * Calculate the equation of center for the sun. This value is a correction
     * to add to the geometric mean longitude in order to get the "true" longitude
     * of the sun.
     *
     * @param  t number of Julian centuries since J2000.
     * @return Equation of center in degrees.
     */
    private static double sunEquationOfCenter(final double t) {
        final double m = Math.toRadians(sunGeometricMeanAnomaly(t));
        return Math.sin(1*m) * (1.914602 - t*(0.004817 + 0.000014*t)) +
               Math.sin(2*m) * (0.019993 - t*(0.000101             )) +
               Math.sin(3*m) * (0.000289);
    }

    /**
     * Calculate the Geometric Mean Longitude of the Sun.
     * This value is close to 0° at the spring equinox,
     * 90° at the summer solstice, 180° at the automne equinox
     * and 270° at the winter solstice.
     *
     * @param  t number of Julian centuries since J2000.
     * @return Geometric Mean Longitude of the Sun in degrees,
     *         in the range 0° (inclusive) to 360° (exclusive).
     */
    private static double sunGeometricMeanLongitude(final double t) {
        double L0 = 280.46646 + t*(36000.76983 + 0.0003032*t);
        L0 = L0 - 360*Math.floor(L0/360);
        return L0;
    }

    /**
     * Calculate the true longitude of the sun. This the geometric mean
     * longitude plus a correction factor ("equation of center" for the
     * sun).
     *
     * @param  t number of Julian centuries since J2000.
     * @return Sun's true longitude in degrees.
     */
    private static double sunTrueLongitude(final double t) {
        return sunGeometricMeanLongitude(t) + sunEquationOfCenter(t);
    }

    /**
     * Calculate the apparent longitude of the sun.
     *
     * @param  t number of Julian centuries since J2000.
     * @return Sun's apparent longitude in degrees.
     */
    private static double sunApparentLongitude(final double t) {
        final double omega = Math.toRadians(125.04 - 1934.136 * t);
        return sunTrueLongitude(t) - 0.00569 - 0.00478 * Math.sin(omega);
    }

    /**
     * Calculate the Geometric Mean Anomaly of the Sun.
     *
     * @param  t number of Julian centuries since J2000.
     * @return Geometric Mean Anomaly of the Sun in degrees.
     */
    private static double sunGeometricMeanAnomaly(final double t) {
        return 357.52911 + t * (35999.05029 - 0.0001537*t);
    }

    /**
     * Calculate the true anamoly of the sun.
     *
     * @param  t number of Julian centuries since J2000.
     * @return Sun's true anamoly in degrees.
     */
    private static double sunTrueAnomaly(final double t) {
        return sunGeometricMeanAnomaly(t) + sunEquationOfCenter(t);
    }

    /**
     * Calculate the eccentricity of earth's orbit. This is the ratio
     * {@code (a-b)/a} where <var>a</var> is the semi-major axis
     * length and <var>b</var> is the semi-minor axis length.   Value
     * is 0 for a circular orbit.
     *
     * @param  t number of Julian centuries since J2000.
     * @return The unitless eccentricity.
     */
    private static double eccentricityEarthOrbit(final double t) {
        return 0.016708634 - t*(0.000042037 + 0.0000001267*t);
    }

    /**
     * Calculate the distance to the sun in Astronomical Units (AU).
     *
     * @param  t number of Julian centuries since J2000.
     * @return Sun radius vector in AUs.
     */
    private static double sunRadiusVector(final double t) {
        final double v = Math.toRadians(sunTrueAnomaly(t));
        final double e = eccentricityEarthOrbit(t);
        return (1.000001018 * (1 - e*e)) / (1 + e*Math.cos(v));
    }

    /**
     * Calculate the mean obliquity of the ecliptic.
     *
     * @param  t number of Julian centuries since J2000.
     * @return Mean obliquity in degrees.
     */
    private static double meanObliquityOfEcliptic(final double t) {
        final double seconds = 21.448 - t*(46.8150 + t*(0.00059 - t*(0.001813)));
        return 23.0 + (26.0 + (seconds/60.0))/60.0;
    }


    /**
     * Calculate the corrected obliquity of the ecliptic.
     *
     * @param  t number of Julian centuries since J2000.
     * @return Corrected obliquity in degrees.
     */
    private static double obliquityCorrected(final double t) {
        final double e0 = meanObliquityOfEcliptic(t);
        final double omega = Math.toRadians(125.04 - 1934.136*t);
        return e0 + 0.00256 * Math.cos(omega);
    }

    /**
     * Calculate the right ascension of the sun. Similar to the angular system
     * used to define latitude and longitude on Earth's surface, right ascension
     * is roughly analogous to longitude, and defines an angular offset from the
     * meridian of the vernal equinox.
     *
     * <P align="center"><img src="doc-files/CelestialSphere.png"></P>
     *
     * @param t number of Julian centuries since J2000.
     * @return Sun's right ascension in degrees.
     */
    private static double sunRightAscension(final double t) {
        final double e = Math.toRadians(obliquityCorrected(t));
        final double b = Math.toRadians(sunApparentLongitude(t));
        final double y = Math.sin(b) * Math.cos(e);
        final double x = Math.cos(b);
        final double alpha = Math.atan2(y, x);
        return Math.toDegrees(alpha);
    }

    /**
     * Calculate the declination of the sun. Declination is analogous to latitude
     * on Earth's surface, and measures an angular displacement north or south
     * from the projection of Earth's equator on the celestial sphere to the
     * location of a celestial body.
     *
     * @param t number of Julian centuries since J2000.
     * @return Sun's declination in degrees.
     */
    private static double sunDeclination(final double t) {
        final double e = Math.toRadians(obliquityCorrected(t));
        final double b = Math.toRadians(sunApparentLongitude(t));
        final double sint = Math.sin(e) * Math.sin(b);
        final double theta = Math.asin(sint);
        return Math.toDegrees(theta);
    }

   /**
    * Calculate the Universal Coordinated Time (UTC) of solar noon for the given day
    * at the given location on earth.
    *
    * @param  lon       longitude of observer in degrees.                               
    * @param  eqTime    Equation of time.
    * @return Time in minutes from beginnning of day in UTC.
    */
    private static double solarNoonTime(final double lon, final double eqTime) { 
        return 720.0 + (lon * 4.0) - eqTime;
    }

    /**
     * Calculate the difference between true solar time and mean. The "equation
     * of time" is a term accounting for changes in the time of solar noon for
     * a given location over the course of a year. Earth's elliptical orbit and
     * Kepler's law of equal areas in equal times are the culprits behind this
     * phenomenon. See the
     * <A HREF="http://www.analemma.com/Pages/framesPage.html">Analemma page</A>.
     * Below is a plot of the equation of time versus the day of the year.
     *
     * <P align="center"><img src="doc-files/EquationOfTime.png"></P>
     *
     * @param  t number of Julian centuries since J2000.
     * @return Equation of time in minutes of time.
     */
    private static double equationOfTime(final double t) {
        double eps = Math.toRadians(obliquityCorrected(t));
        double l0  = Math.toRadians(sunGeometricMeanLongitude(t));
        double m   = Math.toRadians(sunGeometricMeanAnomaly(t));
        double e   = eccentricityEarthOrbit(t);
        double y   = Math.tan(eps/2);
        y *= y;

        double sin2l0 = Math.sin(2 * l0);
        double cos2l0 = Math.cos(2 * l0);
        double sin4l0 = Math.sin(4 * l0);
        double sin1m  = Math.sin(m);
        double sin2m  = Math.sin(2 * m);

        double etime = y*sin2l0 - 2*e*sin1m + 4*e*y*sin1m*cos2l0
                       - 0.5*y*y*sin4l0 - 1.25*e*e*sin2m;

        return Math.toDegrees(etime)*4.0;
    }

    /**
     * Computes the refraction correction angle.
     * The effects of the atmosphere vary with atmospheric pressure, humidity
     * and other variables. Therefore the calculation is approximate. Errors
     * can be expected to increase the further away you are from the equator,
     * because the sun rises and sets at a very shallow angle. Small variations
     * in the atmosphere can have a larger effect.
     *
     * @param  zenith The sun zenith angle in degrees.
     * @return The refraction correction in degrees.
     */
    private static double refractionCorrection(final double zenith) {
        final double exoatmElevation = 90 - zenith;
        if (exoatmElevation > 85) {
            return 0;
        }
        final double refractionCorrection; // In minute of degrees
        final double te = Math.tan(Math.toRadians(exoatmElevation));
        if (exoatmElevation > 5.0) {
            refractionCorrection = 58.1/te - 0.07/(te*te*te) + 0.000086/(te*te*te*te*te);
        } else {
            if (exoatmElevation > -0.575) {
                refractionCorrection =  1735.0 + exoatmElevation *
                                       (-518.2 + exoatmElevation *
                                       ( 103.4 + exoatmElevation *
                                       (-12.79 + exoatmElevation *
                                         0.711)));
            } else {
                refractionCorrection = -20.774 / te;
            }
        }
        return refractionCorrection / 3600;
    }

    /**
     * Constructs a sun relative position calculator.
     */
    public SunRelativePosition() {
    }

    /**
     * Constructs a sun relative position calculator with the specified value
     * for the {@linkplain #setTwilight sun elevation at twilight}.
     *
     * @param twilight The new sun elevation at twilight, or {@link Double#NaN}
     *                 if no twilight value should be taken in account.
     * @throws IllegalArgumentException if the twilight value is illegal.
     */
    public SunRelativePosition(final double twilight) throws IllegalArgumentException {
        setTwilight(twilight);
    }

    /**
     * Calculates solar position for the current date, time and location.
     * Results are reported in azimuth and elevation in degrees.
     */
    private void compute() {
        double latitude  = this.latitude;
        double longitude = this.longitude;

        // NOAA convention use positive longitude west, and negative east.
        // Inverse the sign, in order to be closer to OpenGIS convention.
        longitude = -longitude;

        // Compute: 1) Julian day (days ellapsed since January 1, 4723 BC at 12:00 GMT).
        //          2) Time as the centuries ellapsed since January 1, 2000 at 12:00 GMT.
        final double julianDay = Calendar.julianDay(this.time);
        final double time      = (julianDay-2451545)/36525;

        double solarDec = sunDeclination(time);
        double eqTime   = equationOfTime(time);
        this.noonTime   = Math.round(solarNoonTime(longitude, eqTime) * (60*1000)) +
                          (this.time/DAY_MILLIS)*DAY_MILLIS;

        // Formula below use longitude in degrees. Steps are:
        //   1) Extract the time part of the date, in minutes.
        //   2) Apply a correction for longitude and equation of time.
        //   3) Clamp in a 24 hours range (24 hours == 1440 minutes).
        double trueSolarTime = ((julianDay+0.5) - Math.floor(julianDay+0.5)) * 1440;
        trueSolarTime += (eqTime - 4.0*longitude); // Correction in minutes.
        trueSolarTime -= 1440*Math.floor(trueSolarTime/1440);

        // Convert all angles to radians.  From this point until
        // the end of this method, local variables are always in
        // radians. Output variables ('azimuth' and 'elevation')
        // will still computed in degrees.
        longitude = Math.toRadians(longitude);
        latitude  = Math.toRadians(latitude );
        solarDec  = Math.toRadians(solarDec );

        double csz = Math.sin(latitude) *
                     Math.sin(solarDec) +
                     Math.cos(latitude) *
                     Math.cos(solarDec) *
                     Math.cos(Math.toRadians(trueSolarTime/4 - 180));
        if (csz > +1) csz = +1;
        if (csz < -1) csz = -1;

        final double zenith  = Math.acos(csz);
        final double azDenom = Math.cos(latitude) * Math.sin(zenith);

        //////////////////////////////////////////
        ////    Compute azimuth in degrees    ////
        //////////////////////////////////////////
        if (Math.abs(azDenom) > 0.001) {
            double azRad = ((Math.sin(latitude)*Math.cos(zenith)) - Math.sin(solarDec)) / azDenom;
            if (azRad > +1) azRad = +1;
            if (azRad < -1) azRad = -1;

            azimuth = 180 - Math.toDegrees(Math.acos(azRad));
            if (trueSolarTime > 720) { // 720 minutes == 12 hours
                azimuth = -azimuth;
            }
        } else {
            azimuth = (latitude>0) ? 180 : 0;
        }
        azimuth -= 360*Math.floor(azimuth/360);

        ////////////////////////////////////////////
        ////    Compute elevation in degrees    ////
        ////////////////////////////////////////////
        final double refractionCorrection = refractionCorrection(Math.toDegrees(zenith));
        final double solarZen = Math.toDegrees(zenith) - refractionCorrection;

        elevation = 90 - solarZen;
        if (elevation < twilight) {
            // do not report azimuth & elevation after twilight
            azimuth   = DARK;
            elevation = DARK;
        }
        updated = true;
    }

    /**
     * Set the geographic coordinate where to compute the {@linkplain #getElevation elevation}
     * and {@linkplain #getAzimuth azimuth}.
     *
     * @param longitude The longitude in degrees. Positive values are East; negative values are West.
     * @param latitude  The latitude in degrees. Positive values are North, negative values are South.
     */
    public void setCoordinate(double longitude, double latitude) {
        if (latitude > +89.8) latitude = +89.8;
        if (latitude < -89.8) latitude = -89.8;
        if (latitude != this.latitude || longitude != this.longitude) {
            this.latitude  = latitude;
            this.longitude = longitude;
            this.updated   = false;
        }
    }

    /**
     * Set the geographic coordinate where to compute the {@linkplain #getElevation elevation}
     * and {@linkplain #getAzimuth azimuth}.
     *
     * @param point The geographic coordinates in degrees of longitude and latitude.
     */
    public void setCoordinate(final Point2D point) {
        setCoordinate(point.getX(), point.getY());
    }

    /**
     * Returns the coordinate used for {@linkplain #getElevation elevation} and
     * {@linkplain #getAzimuth azimuth} computation. This is the coordinate
     * specified during the last call to a {@link #setCoordinate(double,double)
     * setCoordinate(...)} method.
     */
    public Point2D getCoordinate() {
        return new Point2D.Double(longitude, latitude);
    }

    /**
     * Set the date and time when to compute the {@linkplain #getElevation elevation}
     * and {@linkplain #getAzimuth azimuth}.
     *
     * @param date The date and time.
     */
    public void setDate(final Date date) {
        final long time = date.getTime();
        if (time != this.time) {
            this.time = time;
            this.updated = false;
        }
    }

    /**
     * Returns the date used for {@linkplain #getElevation elevation} and
     * {@linkplain #getAzimuth azimuth} computation. This is the date specified
     * during the last call to {@link #setDate}.
     */
    public Date getDate() {
        return new Date(time);
    }

    /**
     * Set the sun's {@linkplain #getElevation elevation angle} at twilight, in degrees.
     * Common values are defined for the
     * {@linkplain #ASTRONOMICAL_TWILIGHT astronomical twilight} (-18°),
     * {@linkplain #NAUTICAL_TWILIGHT nautical twilight} (-12°) and
     * {@linkplain #CIVIL_TWILIGHT civil twilight} (-6°).
     * The {@linkplain #getElevation elevation} and {@linkplain #getAzimuth azimuth} are
     * set to {@linkplain Double#NaN NaN} when the sun elevation is below the twilight
     * value (i.e. during night). The default value is {@link #CIVIL_TWILIGHT}.
     *
     * @param twilight The new sun elevation at twilight, or {@link Double#NaN}
     *                 if no twilight value should be taken in account.
     * @throws IllegalArgumentException if the twilight value is illegal.
     */
    public void setTwilight(final double twilight) throws IllegalArgumentException {
        if (twilight<-90 || twilight>-90) {
            // TODO: provides a better (localized) message.
            throw new IllegalArgumentException(String.valueOf(twilight));
        }
        this.twilight = twilight;
        this.updated = false;
    }

    /**
     * Returns the sun's {@linkplain #getElevation elevation angle} at twilight, in degrees.
     * This is the value set during the last call to {@link #setTwilight}.
     */
    public double getTwilight() {
        return twilight;
    }

    /**
     * Retourne l'azimuth en degrés.
     *
     * @return L'azimuth en degrés.
     */
    public double getAzimuth() {
        if (!updated) {
            compute();
        }
        return azimuth;
    }

    /**
     * Retourne l'élévation en degrés.
     *
     * @return L'élévation en degrés.
     */
    public double getElevation() {
        if (!updated) {
            compute();
        }
        return elevation;
    }
    
    /**
     * 
     * @return solar zenith angle in degrees 
     */
    public double getSolarZenith() {
       if (!updated) {
          compute();
       }
       return 90.0 - elevation;
    }

    /**
     * Retourne l'heure à laquelle le soleil est au plus haut. L'heure est
     * retournée en nombre de millisecondes écoulées depuis le debut de la
     * journée (minuit) en heure UTC.
     */
    public long getNoonTime() {
        if (!updated) {
            compute();
        }
        return noonTime % DAY_MILLIS;
    }

    /**
     * Retourne la date à laquelle le soleil est au plus haut dans la journée.
     * Cette méthode est équivalente à {@link #getNoonTime} mais inclue le jour
     * de la date qui avait été spécifiée à la méthode {@link #compute}.
     */
    public Date getNoonDate() {
        if (!updated) {
            compute();
        }
        return new Date(noonTime);
    }

    /**
     * Affiche la position du soleil à la date et coordonnées spécifiée.
     * Cette application peut être lancée avec la syntaxe suivante:
     *
     * <pre>SunRelativePosition <var>[longitude]</var> <var>[latitude]</var> <var>[date]</var></pre>
     *
     * où <var>date</var> est un argument optionel spécifiant la date et l'heure.
     * Si cet argument est omis, la date et heure actuelles seront utilisées.
     */
    public static void main(final String[] args) throws ParseException {
        final DateFormat format=DateFormat.getDateTimeInstance(DateFormat.SHORT, DateFormat.SHORT);
        format.setTimeZone(TimeZone.getTimeZone("UTC"));
        double longitude = 0;
        double latitude  = 0;
        Date time = new Date();
        switch (args.length) {
            case 3: time      = format.parse      (args[2]); // fall through
            case 2: latitude  = Double.parseDouble(args[1]); // fall through
            case 1: longitude = Double.parseDouble(args[0]); // fall through
        }
        final SunRelativePosition calculator = new SunRelativePosition();
        calculator.setDate(time);
        calculator.setCoordinate(longitude, latitude);
        System.out.print("Date (UTC): "); System.out.println(format.format(time));
        System.out.print("Longitude:  "); System.out.println(longitude);
        System.out.print("Latitude:   "); System.out.println(latitude);
        System.out.print("Elevation:  "); System.out.println(calculator.getElevation());
        System.out.print("Azimuth:    "); System.out.println(calculator.getAzimuth());
        System.out.print("Noon date:  "); System.out.println(format.format(calculator.getNoonDate()));
    }
}