001/* 002 * Zmanim Java API 003 * Copyright (C) 2004-2011 Eliyahu Hershfeld 004 * 005 * This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General 006 * Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) 007 * any later version. 008 * 009 * This library is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY; without even the implied 010 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more 011 * details. 012 * You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to 013 * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA, 014 * or connect to: http://www.gnu.org/licenses/old-licenses/lgpl-2.1.html 015 */ 016package net.sourceforge.zmanim.util; 017 018import java.util.Calendar; 019 020/** 021 * Implementation of sunrise and sunset methods to calculate astronomical times. This calculator uses the Java algorithm 022 * written by <a href="http://web.archive.org/web/20090531215353/http://www.kevinboone.com/suntimes.html">Kevin 023 * Boone</a> that is based on the <a href = "http://aa.usno.navy.mil/">US Naval Observatory's</a><a 024 * href="http://aa.usno.navy.mil/publications/docs/asa.php">Almanac</a> for Computer algorithm ( <a 025 * href="http://www.amazon.com/exec/obidos/tg/detail/-/0160515106/">Amazon</a>, <a 026 * href="http://search.barnesandnoble.com/booksearch/isbnInquiry.asp?isbn=0160515106">Barnes & Noble</a>) and is 027 * used with his permission. Added to Kevin's code is adjustment of the zenith to account for elevation. 028 * 029 * @author © Eliyahu Hershfeld 2004 - 2011 030 * @author © Kevin Boone 2000 031 * @version 1.1 032 */ 033public class SunTimesCalculator extends AstronomicalCalculator { 034 /** 035 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getCalculatorName() 036 */ 037 public String getCalculatorName() { 038 return "US Naval Almanac Algorithm"; 039 } 040 041 /** 042 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunrise(Calendar, GeoLocation, double, boolean) 043 */ 044 public double getUTCSunrise(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) { 045 double doubleTime = Double.NaN; 046 047 double elevation = adjustForElevation ? geoLocation.getElevation() : 0; 048 double adjustedZenith = adjustZenith(zenith, elevation); 049 050 doubleTime = getTimeUTC(calendar.get(Calendar.YEAR), calendar.get(Calendar.MONTH) + 1, 051 calendar.get(Calendar.DAY_OF_MONTH), geoLocation.getLongitude(), geoLocation.getLatitude(), 052 adjustedZenith, true); 053 return doubleTime; 054 } 055 056 /** 057 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunset(Calendar, GeoLocation, double, boolean) 058 */ 059 public double getUTCSunset(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) { 060 double doubleTime = Double.NaN; 061 double elevation = adjustForElevation ? geoLocation.getElevation() : 0; 062 double adjustedZenith = adjustZenith(zenith, elevation); 063 064 doubleTime = getTimeUTC(calendar.get(Calendar.YEAR), calendar.get(Calendar.MONTH) + 1, 065 calendar.get(Calendar.DAY_OF_MONTH), geoLocation.getLongitude(), geoLocation.getLatitude(), 066 adjustedZenith, false); 067 return doubleTime; 068 } 069 070 /** 071 * The number of degrees of longitude that corresponds to one hour time difference. 072 */ 073 private static final double DEG_PER_HOUR = 360.0 / 24.0; 074 075 /** 076 * sin of an angle in degrees 077 */ 078 private static double sinDeg(double deg) { 079 return Math.sin(deg * 2.0 * Math.PI / 360.0); 080 } 081 082 /** 083 * acos of an angle, result in degrees 084 */ 085 private static double acosDeg(double x) { 086 return Math.acos(x) * 360.0 / (2 * Math.PI); 087 } 088 089 /** 090 * asin of an angle, result in degrees 091 */ 092 private static double asinDeg(double x) { 093 return Math.asin(x) * 360.0 / (2 * Math.PI); 094 } 095 096 /** 097 * tan of an angle in degrees 098 */ 099 private static double tanDeg(double deg) { 100 return Math.tan(deg * 2.0 * Math.PI / 360.0); 101 } 102 103 /** 104 * cos of an angle in degrees 105 */ 106 private static double cosDeg(double deg) { 107 return Math.cos(deg * 2.0 * Math.PI / 360.0); 108 } 109 110 /** 111 * Calculate the day of the year, where Jan 1st is day 1. Note that this method needs to know the year, because leap 112 * years have an impact here. Returns identical output to {@code Calendar.get(Calendar.DAY_OF_YEAR)}. 113 */ 114 private static int getDayOfYear(int year, int month, int day) { 115 int n1 = 275 * month / 9; 116 int n2 = (month + 9) / 12; 117 int n3 = (1 + ((year - 4 * (year / 4) + 2) / 3)); 118 return n1 - (n2 * n3) + day - 30; 119 } 120 121 /** 122 * Get time difference between location's longitude and the Meridian, in hours. West of Meridian has a negative time 123 * difference 124 */ 125 private static double getHoursFromMeridian(double longitude) { 126 return longitude / DEG_PER_HOUR; 127 } 128 129 /** 130 * Gets the approximate time of sunset or sunrise In _days_ since midnight Jan 1st, assuming 6am and 6pm events. We 131 * need this figure to derive the Sun's mean anomaly 132 */ 133 private static double getApproxTimeDays(int dayOfYear, double hoursFromMeridian, boolean isSunrise) { 134 if (isSunrise) { 135 return dayOfYear + ((6.0 - hoursFromMeridian) / 24); 136 } else { // sunset 137 return dayOfYear + ((18.0 - hoursFromMeridian) / 24); 138 } 139 } 140 141 /** 142 * Calculate the Sun's mean anomaly in degrees, at sunrise or sunset, given the longitude in degrees 143 */ 144 private static double getMeanAnomaly(int dayOfYear, double longitude, boolean isSunrise) { 145 return (0.9856 * getApproxTimeDays(dayOfYear, getHoursFromMeridian(longitude), isSunrise)) - 3.289; 146 } 147 148 /** 149 * Calculates the Sun's true longitude in degrees. The result is an angle gte 0 and lt 360. Requires the Sun's mean 150 * anomaly, also in degrees 151 */ 152 private static double getSunTrueLongitude(double sunMeanAnomaly) { 153 double l = sunMeanAnomaly + (1.916 * sinDeg(sunMeanAnomaly)) + (0.020 * sinDeg(2 * sunMeanAnomaly)) + 282.634; 154 155 // get longitude into 0-360 degree range 156 if (l >= 360.0) { 157 l = l - 360.0; 158 } 159 if (l < 0) { 160 l = l + 360.0; 161 } 162 return l; 163 } 164 165 /** 166 * Calculates the Sun's right ascension in hours, given the Sun's true longitude in degrees. Input and output are 167 * angles gte 0 and lt 360. 168 */ 169 private static double getSunRightAscensionHours(double sunTrueLongitude) { 170 double a = 0.91764 * tanDeg(sunTrueLongitude); 171 double ra = 360.0 / (2.0 * Math.PI) * Math.atan(a); 172 173 double lQuadrant = Math.floor(sunTrueLongitude / 90.0) * 90.0; 174 double raQuadrant = Math.floor(ra / 90.0) * 90.0; 175 ra = ra + (lQuadrant - raQuadrant); 176 177 return ra / DEG_PER_HOUR; // convert to hours 178 } 179 180 /** 181 * Gets the cosine of the Sun's local hour angle 182 */ 183 private static double getCosLocalHourAngle(double sunTrueLongitude, double latitude, double zenith) { 184 double sinDec = 0.39782 * sinDeg(sunTrueLongitude); 185 double cosDec = cosDeg(asinDeg(sinDec)); 186 return (cosDeg(zenith) - (sinDec * sinDeg(latitude))) / (cosDec * cosDeg(latitude)); 187 } 188 189 /** 190 * Calculate local mean time of rising or setting. By `local' is meant the exact time at the location, assuming that 191 * there were no time zone. That is, the time difference between the location and the Meridian depended entirely on 192 * the longitude. We can't do anything with this time directly; we must convert it to UTC and then to a local time. 193 * The result is expressed as a fractional number of hours since midnight 194 */ 195 private static double getLocalMeanTime(double localHour, double sunRightAscensionHours, double approxTimeDays) { 196 return localHour + sunRightAscensionHours - (0.06571 * approxTimeDays) - 6.622; 197 } 198 199 /** 200 * Get sunrise or sunset time in UTC, according to flag. 201 * 202 * @param year 203 * 4-digit year 204 * @param month 205 * month, 1-12 (not the zero based Java month 206 * @param day 207 * day of month, 1-31 208 * @param longitude 209 * in degrees, longitudes west of Meridian are negative 210 * @param latitude 211 * in degrees, latitudes south of equator are negative 212 * @param zenith 213 * Sun's zenith, in degrees 214 * @param type 215 * type of calculation to carry out {@link #TYPE_SUNRISE} or {@link #TYPE_SUNRISE}. 216 * 217 * @return the time as a double. If an error was encountered in the calculation (expected behavior for some 218 * locations such as near the poles, {@link Double.NaN} will be returned. 219 */ 220 private static double getTimeUTC(int year, int month, int day, double longitude, double latitude, double zenith, 221 boolean isSunrise) { 222 int dayOfYear = getDayOfYear(year, month, day); 223 double sunMeanAnomaly = getMeanAnomaly(dayOfYear, longitude, isSunrise); 224 double sunTrueLong = getSunTrueLongitude(sunMeanAnomaly); 225 double sunRightAscensionHours = getSunRightAscensionHours(sunTrueLong); 226 double cosLocalHourAngle = getCosLocalHourAngle(sunTrueLong, latitude, zenith); 227 228 double localHourAngle = 0; 229 if (isSunrise) { 230 localHourAngle = 360.0 - acosDeg(cosLocalHourAngle); 231 } else { // sunset 232 localHourAngle = acosDeg(cosLocalHourAngle); 233 } 234 double localHour = localHourAngle / DEG_PER_HOUR; 235 236 double localMeanTime = getLocalMeanTime(localHour, sunRightAscensionHours, 237 getApproxTimeDays(dayOfYear, getHoursFromMeridian(longitude), isSunrise)); 238 double pocessedTime = localMeanTime - getHoursFromMeridian(longitude); 239 while (pocessedTime < 0.0) { 240 pocessedTime += 24.0; 241 } 242 while (pocessedTime >= 24.0) { 243 pocessedTime -= 24.0; 244 } 245 return pocessedTime; 246 } 247}