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