// 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));
}
}