001 /* 002 * Zmanim Java API 003 * Copyright (C) 2004-2011 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 - 2011 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 this.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 adjustedZenith = zenith; 053 double doubleTime = Double.NaN; 054 055 if (adjustForElevation) { 056 adjustedZenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation() 057 .getElevation()); 058 } else { 059 adjustedZenith = adjustZenith(zenith, 0); 060 } 061 doubleTime = getTimeUTC(astronomicalCalendar.getCalendar().get( 062 Calendar.YEAR), astronomicalCalendar.getCalendar().get( 063 Calendar.MONTH) + 1, astronomicalCalendar.getCalendar().get( 064 Calendar.DAY_OF_MONTH), astronomicalCalendar.getGeoLocation() 065 .getLongitude(), astronomicalCalendar.getGeoLocation() 066 .getLatitude(), adjustedZenith, TYPE_SUNRISE); 067 return doubleTime; 068 } 069 070 /** 071 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunset(AstronomicalCalendar, 072 * double, boolean) 073 */ 074 public double getUTCSunset(AstronomicalCalendar astronomicalCalendar, 075 double zenith, boolean adjustForElevation) { 076 double doubleTime = Double.NaN; 077 double adjustedZenith = zenith; 078 079 if (adjustForElevation) { 080 adjustedZenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation() 081 .getElevation()); 082 } else { 083 adjustedZenith = adjustZenith(zenith, 0); 084 } 085 doubleTime = getTimeUTC(astronomicalCalendar.getCalendar().get( 086 Calendar.YEAR), astronomicalCalendar.getCalendar().get( 087 Calendar.MONTH) + 1, astronomicalCalendar.getCalendar().get( 088 Calendar.DAY_OF_MONTH), astronomicalCalendar.getGeoLocation() 089 .getLongitude(), astronomicalCalendar.getGeoLocation() 090 .getLatitude(), adjustedZenith, TYPE_SUNSET); 091 return doubleTime; 092 } 093 094 /** 095 * Default value for Sun's zenith and true rise/set 096 */ 097 public static final double ZENITH = 90 + 50.0 / 60.0; 098 099 private static final int TYPE_SUNRISE = 0; 100 101 private static final int TYPE_SUNSET = 1; 102 103 // DEG_PER_HOUR is the number of degrees of longitude 104 // that corresponds to one hour time difference. 105 private static final double DEG_PER_HOUR = 360.0 / 24.0; 106 107 /** 108 * sin of an angle in degrees 109 */ 110 private static double sinDeg(double deg) { 111 return Math.sin(deg * 2.0 * Math.PI / 360.0); 112 } 113 114 /** 115 * acos of an angle, result in degrees 116 */ 117 private static double acosDeg(double x) { 118 return Math.acos(x) * 360.0 / (2 * Math.PI); 119 } 120 121 /** 122 * asin of an angle, result in degrees 123 */ 124 private static double asinDeg(double x) { 125 return Math.asin(x) * 360.0 / (2 * Math.PI); 126 } 127 128 /** 129 * tan of an angle in degrees 130 */ 131 private static double tanDeg(double deg) { 132 return Math.tan(deg * 2.0 * Math.PI / 360.0); 133 } 134 135 /** 136 * cos of an angle in degrees 137 */ 138 private static double cosDeg(double deg) { 139 return Math.cos(deg * 2.0 * Math.PI / 360.0); 140 } 141 142 /** 143 * Calculate the day of the year, where Jan 1st is day 1. Note that this 144 * method needs to know the year, because leap years have an impact here 145 */ 146 private static int getDayOfYear(int year, int month, int day) { 147 int n1 = 275 * month / 9; 148 int n2 = (month + 9) / 12; 149 int n3 = (1 + ((year - 4 * (year / 4) + 2) / 3)); 150 int n = n1 - (n2 * n3) + day - 30; 151 return n; 152 } 153 154 /** 155 * Get time difference between location's longitude and the Meridian, in 156 * hours. West of Meridian has a negative time difference 157 */ 158 private static double getHoursFromMeridian(double longitude) { 159 return longitude / DEG_PER_HOUR; 160 } 161 162 /** 163 * Gets the approximate time of sunset or sunrise In _days_ since midnight 164 * Jan 1st, assuming 6am and 6pm events. We need this figure to derive the 165 * Sun's mean anomaly 166 */ 167 private static double getApproxTimeDays(int dayOfYear, 168 double hoursFromMeridian, int type) { 169 if (type == TYPE_SUNRISE) { 170 return dayOfYear + ((6.0 - hoursFromMeridian) / 24); 171 } else /* if (type == TYPE_SUNSET) */{ 172 return dayOfYear + ((18.0 - hoursFromMeridian) / 24); 173 } 174 } 175 176 /** 177 * Calculate the Sun's mean anomaly in degrees, at sunrise or sunset, given 178 * the longitude in degrees 179 */ 180 private static double getMeanAnomaly(int dayOfYear, double longitude, 181 int type) { 182 return (0.9856 * getApproxTimeDays(dayOfYear, 183 getHoursFromMeridian(longitude), type)) - 3.289; 184 } 185 186 /** 187 * Calculates the Sun's true longitude in degrees. The result is an angle 188 * gte 0 and lt 360. Requires the Sun's mean anomaly, also in degrees 189 */ 190 private static double getSunTrueLongitude(double sunMeanAnomaly) { 191 double l = sunMeanAnomaly + (1.916 * sinDeg(sunMeanAnomaly)) 192 + (0.020 * sinDeg(2 * sunMeanAnomaly)) + 282.634; 193 194 // get longitude into 0-360 degree range 195 if (l >= 360.0) { 196 l = l - 360.0; 197 } 198 if (l < 0) { 199 l = l + 360.0; 200 } 201 return l; 202 } 203 204 /** 205 * Calculates the Sun's right ascension in hours, given the Sun's true 206 * longitude in degrees. Input and output are angles gte 0 and lt 360. 207 */ 208 private static double getSunRightAscensionHours(double sunTrueLongitude) { 209 double a = 0.91764 * tanDeg(sunTrueLongitude); 210 double ra = 360.0 / (2.0 * Math.PI) * Math.atan(a); 211 // get result into 0-360 degree range 212 // if (ra >= 360.0) ra = ra - 360.0; 213 // if (ra < 0) ra = ra + 360.0; 214 215 double lQuadrant = Math.floor(sunTrueLongitude / 90.0) * 90.0; 216 double raQuadrant = Math.floor(ra / 90.0) * 90.0; 217 ra = ra + (lQuadrant - raQuadrant); 218 219 return ra / DEG_PER_HOUR; // convert to hours 220 } 221 222 /** 223 * Gets the cosine of the Sun's local hour angle 224 */ 225 private static double getCosLocalHourAngle(double sunTrueLongitude, 226 double latitude, double zenith) { 227 double sinDec = 0.39782 * sinDeg(sunTrueLongitude); 228 double cosDec = cosDeg(asinDeg(sinDec)); 229 230 double cosH = (cosDeg(zenith) - (sinDec * sinDeg(latitude))) 231 / (cosDec * cosDeg(latitude)); 232 233 // Check bounds 234 235 return cosH; 236 } 237 238 /** 239 * Gets the cosine of the Sun's local hour angle for default zenith 240 */ 241 // private static double getCosLocalHourAngle(double sunTrueLongitude, 242 // double latitude) { 243 // return getCosLocalHourAngle(sunTrueLongitude, latitude, ZENITH); 244 // } 245 246 /** 247 * Calculate local mean time of rising or setting. By `local' is meant the 248 * exact time at the location, assuming that there were no time zone. That 249 * is, the time difference between the location and the Meridian depended 250 * entirely on the longitude. We can't do anything with this time directly; 251 * we must convert it to UTC and then to a local time. The result is 252 * expressed as a fractional number of hours since midnight 253 */ 254 private static double getLocalMeanTime(double localHour, 255 double sunRightAscensionHours, double approxTimeDays) { 256 return localHour + sunRightAscensionHours - (0.06571 * approxTimeDays) 257 - 6.622; 258 } 259 260 /** 261 * Get sunrise or sunset time in UTC, according to flag. 262 * 263 * @param year 264 * 4-digit year 265 * @param month 266 * month, 1-12 (not the zero based Java month 267 * @param day 268 * day of month, 1-31 269 * @param longitude 270 * in degrees, longitudes west of Meridian are negative 271 * @param latitude 272 * in degrees, latitudes south of equator are negative 273 * @param zenith 274 * Sun's zenith, in degrees 275 * @param type 276 * type of calculation to carry out {@link #TYPE_SUNRISE} or 277 * {@link #TYPE_SUNRISE}. 278 * 279 * @return the time as a double. If an error was encountered in the 280 * calculation (expected behavior for some locations such as near 281 * the poles, {@link Double.NaN} will be returned. 282 */ 283 private static double getTimeUTC(int year, int month, int day, 284 double longitude, double latitude, double zenith, int type) { 285 int dayOfYear = getDayOfYear(year, month, day); 286 double sunMeanAnomaly = getMeanAnomaly(dayOfYear, longitude, type); 287 double sunTrueLong = getSunTrueLongitude(sunMeanAnomaly); 288 double sunRightAscensionHours = getSunRightAscensionHours(sunTrueLong); 289 double cosLocalHourAngle = getCosLocalHourAngle(sunTrueLong, latitude, 290 zenith); 291 292 double localHourAngle = 0; 293 if (type == TYPE_SUNRISE) { 294 if (cosLocalHourAngle > 1) { // no rise. No need for an Exception 295 // since the calculation 296 // will return Double.NaN 297 } 298 localHourAngle = 360.0 - acosDeg(cosLocalHourAngle); 299 } else /* if (type == TYPE_SUNSET) */{ 300 if (cosLocalHourAngle < -1) {// no SET. No need for an Exception 301 // since the calculation 302 // will return Double.NaN 303 } 304 localHourAngle = acosDeg(cosLocalHourAngle); 305 } 306 double localHour = localHourAngle / DEG_PER_HOUR; 307 308 double localMeanTime = getLocalMeanTime(localHour, 309 sunRightAscensionHours, getApproxTimeDays(dayOfYear, 310 getHoursFromMeridian(longitude), type)); 311 double pocessedTime = localMeanTime - getHoursFromMeridian(longitude); 312 while (pocessedTime < 0.0) { 313 pocessedTime += 24.0; 314 } 315 while (pocessedTime >= 24.0) { 316 pocessedTime -= 24.0; 317 } 318 return pocessedTime; 319 } 320 321 }