Package jsky.plot

Source Code of jsky.plot.SunRiseSet

// Copyright 2003 Association for Universities for Research in
// Astronomy, Inc., Observatory Control System, Gemini Telescopes
// Project.
//
// $Id: SunRiseSet.java,v 1.2 2009/02/20 23:10:11 abrighto Exp $

package jsky.plot;

import java.util.Calendar;
import java.util.Date;
import java.util.TimeZone;
import java.text.SimpleDateFormat;

import jsky.coords.SiteDesc;
import jsky.util.CalendarUtil;


/**
* Utility class for calculating the times for sunrise, sunset and twilight for a given location and date.
* Based on the algorithm found <a href="http://williams.best.vwh.net/sunrise_sunset_algorithm.htm">here</a>.
*
* @version $Revision: 1.2 $
* @author Allan Brighton
*/
public class SunRiseSet {

    private Date _date;
    private Date _sunset;
    private Date _sunrise;
    private Date _civilTwilightStart;
    private Date _civilTwilightEnd;
    private Date _nauticalTwilightStart;
    private Date _nauticalTwilightEnd;
    private Date _astronomicalTwilightStart;
    private Date _astronomicalTwilightEnd;


    /**
     * Calculates the times for sunrise, sunset and twilight for the given date and location.
     *
     * @param date (in) the date of interest
     * @param site (in) describes the observatory location
     */
    public SunRiseSet(Date date, SiteDesc site) {
        _date = date;
        double longitude = site.getLongitude();
        double latitude = site.getLatitude();

        Calendar cal = Calendar.getInstance(site.getTimeZone());
        cal.setTime(date);
        int dayOfYear = cal.get(Calendar.DAY_OF_YEAR); // want tonight and tomorrow morning

        // convert the longitude to hour value and calculate an approximate time
        double lngHour = longitude / 15.;

        // calculate the times for nautical and astronomical
        _sunset = _calcTime(false, 90. + 50. / 60., dayOfYear, lngHour, latitude);
        _sunrise = _calcTime(true, 90. + 50. / 60., dayOfYear, lngHour, latitude);
        _civilTwilightStart = _calcTime(false, 96., dayOfYear, lngHour, latitude);
        _civilTwilightEnd = _calcTime(true, 96., dayOfYear, lngHour, latitude);
        _nauticalTwilightStart = _calcTime(false, 102., dayOfYear, lngHour, latitude);
        _nauticalTwilightEnd = _calcTime(true, 102., dayOfYear, lngHour, latitude);
        _astronomicalTwilightStart = _calcTime(false, 108., dayOfYear, lngHour, latitude);
        _astronomicalTwilightEnd = _calcTime(true, 108., dayOfYear, lngHour, latitude);
    }

    public Date getSunset() {
        return _sunset;
    }

    public Date getSunrise() {
        return _sunrise;
    }

    public Date getCivilTwilightStart() {
        return _civilTwilightStart;
    }

    public Date getCivilTwilightEnd() {
        return _civilTwilightEnd;
    }

    public Date getNauticalTwilightStart() {
        return _nauticalTwilightStart;
    }

    public Date getNauticalTwilightEnd() {
        return _nauticalTwilightEnd;
    }

    public Date getAstronomicalTwilightStart() {
        return _astronomicalTwilightStart;
    }

    public Date getAstronomicalTwilightEnd() {
        return _astronomicalTwilightEnd;
    }

    // Shortcut for math functions in degrees
    private static double _sin(double d) {
        return Math.sin(Math.toRadians(d));
    }

    private static double _cos(double d) {
        return Math.cos(Math.toRadians(d));
    }

    private static double _tan(double d) {
        return Math.tan(Math.toRadians(d));
    }

    private static double _asin(double d) {
        return Math.toDegrees(Math.asin(d));
    }

    private static double _acos(double d) {
        return Math.toDegrees(Math.acos(d));
    }

    private static double _atan(double d) {
        return Math.toDegrees(Math.atan(d));
    }


    // Returns the time for sunrise (if rise is true) or sunset (if false),
    // using the given zenith (in degrees):
    //    offical      = 90 degrees 50'
    //    civil        = 96 degrees
    //    nautical     = 102 degrees
    //    astronomical = 108 degrees
    // Assumes sunrise is on next day.
    private Date _calcTime(boolean rise, double zenith, int dayOfYear,
                           double lngHour, double latitude) {

        if (rise) {
            dayOfYear++;
        }
        double t = (rise ? (dayOfYear + ((6. - lngHour) / 24.)) : (dayOfYear + ((18. - lngHour) / 24.)));

        // calculate the Sun's mean anomaly
        double m = (0.9856 * t) - 3.289;

        // calculate the Sun's true longitude
        double l = m + (1.916 * _sin(m)) + (0.020 * _sin(2 * m)) + 282.634;

        // l potentially needs to be adjusted into the range [0,360)
        l = _normalize(l, 360.);

        // calculate the Sun's right ascension
        double ra = _atan(0.91764 * _tan(l));

        // ra potentially needs to be adjusted into the range [0,360)
        ra = _normalize(ra, 360.);

        // right ascension value needs to be in the same quadrant as l
        double lQuadrant = (Math.floor(l / 90)) * 90;
        double raQuadrant = (Math.floor(ra / 90)) * 90;
        ra += (lQuadrant - raQuadrant);

        // right ascension value needs to be converted into hours
        ra /= 15;

        // calculate the Sun's declination
        double sinDec = 0.39782 * _sin(l);
        double cosDec = _cos(_asin(sinDec));

        // calculate the Sun's local hour angle
        double cosH = (_cos(zenith) - (sinDec * _sin(latitude))) / (cosDec * _cos(latitude));
        if (cosH > 1) {
            // the sun never rises on this location (on the specified date)
        } else if (cosH < -1) {
            // the sun never sets on this location (on the specified date)
        }

        // finish calculating H and convert into hours
        double h = (rise ? (360 - _acos(cosH)) : _acos(cosH)) / 15.;

        // calculate local mean time of rising/setting
        t = h + ra - (0.06571 * t) - 6.622;

        // adjust back to UTC
        double ut = t - lngHour;

        // UT potentially needs to be adjusted into the range [0,24) by adding/subtracting 24
//        ut = _normalize(ut, 24.);

        return _getDate(ut);
    }


    // Return a value between 0 and maxValue by possibly adding or subtracting maxValue
    private static double _normalize(double value, double maxValue) {
        if (value > maxValue)
            value -= maxValue;
        else if (value < 0)
            value += maxValue;
        return value;
    }

    // print out the date in the given time zone
    private static String _fmt(Date date, TimeZone tz) {
        SimpleDateFormat dateFormat = new SimpleDateFormat("EEE, d MMM yyyy HH:mm:ss z");
        dateFormat.setTimeZone(tz);
        return dateFormat.format(date);
    }

    // Return a Date object for the given UT hours (0..24).
    // If rise is true, use the next day.
    private Date _getDate(double hours) {
        Calendar cal = Calendar.getInstance(TimeZone.getTimeZone("UT"));
        cal.setTime(_date);
        int h = cal.get(Calendar.HOUR_OF_DAY);
        boolean nextDay = (hours < h);
        CalendarUtil.setHours(cal, hours, nextDay);
        return cal.getTime();
    }


    /** Test main. */
    public static void main(String[] args) {
        SiteDesc site = ElevationPlotUtil.MAUNA_KEA;
        TimeZone tz = site.getTimeZone();
        Calendar cal = Calendar.getInstance(tz);

        cal.set(2005, Calendar.MARCH, 8, 12, 0, 0);
        Date date = cal.getTime();
        SunRiseSet srs = new SunRiseSet(date, site);
        System.out.println("Date = " + _fmt(date, tz) + ", site = " + site.getName());
        System.out.println("Sunset                      = " + _fmt(srs.getSunset(), tz));
        System.out.println("Sunrise                     = " + _fmt(srs.getSunrise(), tz));
        System.out.println("CivilTwilight Start         = " + _fmt(srs.getCivilTwilightStart(), tz));
        System.out.println("Civil Twilight End          = " + _fmt(srs.getCivilTwilightEnd(), tz));
        System.out.println("NauticalTwilight Start      = " + _fmt(srs.getNauticalTwilightStart(), tz));
        System.out.println("NauticalTwilight End        = " + _fmt(srs.getNauticalTwilightEnd(), tz));
        System.out.println("Astronomical Twilight Start = " + _fmt(srs.getAstronomicalTwilightStart(), tz));
        System.out.println("AstronomicalTwilight End    = " + _fmt(srs.getAstronomicalTwilightEnd(), tz));

        System.out.println();
        cal.set(2005, Calendar.MARCH, 9, 12, 0, 0);
        date = cal.getTime();
        srs = new SunRiseSet(date, site);
        tz = site.getTimeZone();
        System.out.println("Date = " + _fmt(date, tz) + ", site = " + site.getName());
        System.out.println("Sunset                      = " + _fmt(srs.getSunset(), tz));
        System.out.println("Sunrise                     = " + _fmt(srs.getSunrise(), tz));
        System.out.println("CivilTwilight Start         = " + _fmt(srs.getCivilTwilightStart(), tz));
        System.out.println("Civil Twilight End          = " + _fmt(srs.getCivilTwilightEnd(), tz));
        System.out.println("NauticalTwilight Start      = " + _fmt(srs.getNauticalTwilightStart(), tz));
        System.out.println("NauticalTwilight End        = " + _fmt(srs.getNauticalTwilightEnd(), tz));
        System.out.println("Astronomical Twilight Start = " + _fmt(srs.getAstronomicalTwilightStart(), tz));
        System.out.println("AstronomicalTwilight End    = " + _fmt(srs.getAstronomicalTwilightEnd(), tz));
    }
}

TOP

Related Classes of jsky.plot.SunRiseSet

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.