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 net.sourceforge.zmanim.AstronomicalCalendar; 020 import java.util.Calendar; 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.jstot.me.uk/jsuntimes/">Jonathan Stott</a> that is based on 026 * the implementation by <a href=""http://noaa.gov">NOAA - National Oceanic and 027 * Atmospheric Administration</a>'s <a href = 028 * "http://www.srrb.noaa.gov/highlights/sunrise/sunrisehtml">Surface Radiation 029 * Research Branch</a>. NOAA's <a 030 * href="http://www.srrb.noaa.gov/highlights/sunrise/solareqns.PDF">implementation</a> 031 * is based on equations from <a 032 * href="http://www.willbell.com/math/mc1.htm">Astronomical Algorithms</a> by 033 * <a href="http://en.wikipedia.org/wiki/Jean_Meeus">Jean Meeus</a>. Jonathan's 034 * implementation was released under the GPL. Added to the algorithm is an 035 * adjustment of the zenith to account for elevation. 036 * 037 * @author Jonathan Stott 2000 - 2004 038 * @author © Eliyahu Hershfeld 2004 - 2008 039 * @version 1.1 040 * @deprecated This class is based on the NOAA algorithm but does not return calculations that match the NOAA algorithm JavaScript implementation. The calculations are about 2 minutes off. This call has been replaced by the NOAACalculator class. 041 * @see net.sourceforge.zmanim.util.NOAACalculator 042 */ 043 public class JSuntimeCalculator extends AstronomicalCalculator { 044 private String calculatorName = "US National Oceanic and Atmospheric Administration Algorithm"; 045 /** 046 * @deprecated This class is based on the NOAA algorithm but does not return calculations that match the NOAA algorithm JavaScript implementation. The calculations are about 2 minutes off. This call has been replaced by the NOAACalculator class. 047 * @see net.sourceforge.zmanim.util.NOAACalculator#getCalculatorName() 048 */ 049 public String getCalculatorName() { 050 return calculatorName; 051 } 052 053 /** 054 * @deprecated This class is based on the NOAA algorithm but does not return calculations that match the NOAA algorithm JavaScript implementation. The calculations are about 2 minutes off. This call has been replaced by the NOAACalculator class. 055 * @see net.sourceforge.zmanim.util.NOAACalculator#getUTCSunrise(AstronomicalCalendar, double, boolean) 056 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunrise(AstronomicalCalendar, 057 * double, boolean) 058 * @throws ZmanimException 059 * if the year entered == 2000. This calculator can't properly 060 * deal with the year 2000. It can properly calculate times for 061 * years <> 2000. 062 */ 063 public double getUTCSunrise(AstronomicalCalendar astronomicalCalendar, double zenith, boolean adjustForElevation) { 064 // if (astronomicalCalendar.getCalendar().get(Calendar.YEAR) == 2000) { 065 // throw new ZmanimException( 066 // "JSuntimeCalculator can not calculate times for the year 2000. Please try a date with a different year."); 067 // } 068 069 if (adjustForElevation) { 070 zenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation() 071 .getElevation()); 072 } else { 073 zenith = adjustZenith(zenith, 0); 074 } 075 double timeMins = morningPhenomenon(dateToJulian(astronomicalCalendar 076 .getCalendar()), astronomicalCalendar.getGeoLocation() 077 .getLatitude(), -astronomicalCalendar.getGeoLocation() 078 .getLongitude(), zenith); 079 return timeMins / 60; 080 } 081 082 /** 083 * @deprecated This class is based on the NOAA algorithm but does not return calculations that match the NOAAA algorithm JavaScript implementation. The calculations are about 2 minutes off. This call has been replaced by the NOAACalculator class. 084 * @see net.sourceforge.zmanim.util.NOAACalculator#getUTCSunset(AstronomicalCalendar, double, boolean) 085 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunset(AstronomicalCalendar, 086 * double, boolean) 087 * @throws ZmanimException 088 * if the year entered == 2000. This calculator can't properly 089 * deal with the year 2000. It can properly calculate times for 090 * years <> 2000. 091 */ 092 public double getUTCSunset(AstronomicalCalendar astronomicalCalendar, double zenith, boolean adjustForElevation) { 093 // if (astronomicalCalendar.getCalendar().get(Calendar.YEAR) == 2000) { 094 // throw new ZmanimException( 095 // "JSuntimeCalculator can not calculate times for the year 2000. Please try a date with a different year."); 096 // } 097 098 if (adjustForElevation) { 099 zenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation() 100 .getElevation()); 101 } else { 102 zenith = adjustZenith(zenith, 0); 103 } 104 double timeMins = eveningPhenomenon(dateToJulian(astronomicalCalendar 105 .getCalendar()), astronomicalCalendar.getGeoLocation() 106 .getLatitude(), -astronomicalCalendar.getGeoLocation() 107 .getLongitude(), zenith); 108 return timeMins / 60; 109 } 110 111 /** 112 * Calculate the UTC of a morning phenomenon for the given day at the given 113 * latitude and longitude on Earth 114 * 115 * @param julian 116 * Julian day 117 * @param latitude 118 * latitude of observer in degrees 119 * @param longitude 120 * longitude of observer in degrees 121 * @param zenithDistance 122 * one of Sun.SUNRISE_SUNSET_ZENITH_DISTANCE, 123 * Sun.CIVIL_TWILIGHT_ZENITH_DISTANCE, 124 * Sun.NAUTICAL_TWILIGHT_ZENITH_DISTANCE, 125 * Sun.ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE. 126 * @return time in minutes from zero Z 127 */ 128 private static double morningPhenomenon(double julian, double latitude, 129 double longitude, double zenithDistance) { 130 double t = julianDayToJulianCenturies(julian); 131 double eqtime = equationOfTime(t); 132 double solarDec = sunDeclination(t); 133 double hourangle = hourAngleMorning(latitude, solarDec, zenithDistance); 134 double delta = longitude - Math.toDegrees(hourangle); 135 double timeDiff = 4 * delta; 136 double timeUTC = 720 + timeDiff - eqtime; 137 138 // Second pass includes fractional julian day in gamma calc 139 double newt = julianDayToJulianCenturies(julianCenturiesToJulianDay(t) 140 + timeUTC / 1440); 141 eqtime = equationOfTime(newt); 142 solarDec = sunDeclination(newt); 143 hourangle = hourAngleMorning(latitude, solarDec, zenithDistance); 144 delta = longitude - Math.toDegrees(hourangle); 145 timeDiff = 4 * delta; 146 147 double morning = 720 + timeDiff - eqtime; 148 return morning; 149 } 150 151 /** 152 * Calculate the UTC of an evening phenomenon for the given day at the given 153 * latitude and longitude on Earth 154 * 155 * @param julian 156 * Julian day 157 * @param latitude 158 * latitude of observer in degrees 159 * @param longitude 160 * longitude of observer in degrees 161 * @param zenithDistance 162 * one of Sun.SUNRISE_SUNSET_ZENITH_DISTANCE, 163 * Sun.CIVIL_TWILIGHT_ZENITH_DISTANCE, 164 * Sun.NAUTICAL_TWILIGHT_ZENITH_DISTANCE, 165 * Sun.ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE. 166 * @return time in minutes from zero Z 167 */ 168 private static double eveningPhenomenon(double julian, double latitude, 169 double longitude, double zenithDistance) { 170 double t = julianDayToJulianCenturies(julian); 171 172 // First calculates sunrise and approx length of day 173 double eqtime = equationOfTime(t); 174 double solarDec = sunDeclination(t); 175 double hourangle = hourAngleEvening(latitude, solarDec, zenithDistance); 176 177 double delta = longitude - Math.toDegrees(hourangle); 178 double timeDiff = 4 * delta; 179 double timeUTC = 720 + timeDiff - eqtime; 180 181 // first pass used to include fractional day in gamma calc 182 double newt = julianDayToJulianCenturies(julianCenturiesToJulianDay(t) 183 + timeUTC / 1440); 184 eqtime = equationOfTime(newt); 185 solarDec = sunDeclination(newt); 186 hourangle = hourAngleEvening(latitude, solarDec, zenithDistance); 187 188 delta = longitude - Math.toDegrees(hourangle); 189 timeDiff = 4 * delta; 190 191 double evening = 720 + timeDiff - eqtime; 192 return evening; 193 } 194 195 private static double dateToJulian(Calendar date) { 196 int year = date.get(Calendar.YEAR); 197 int month = date.get(Calendar.MONTH) + 1; 198 int day = date.get(Calendar.DAY_OF_MONTH); 199 int hour = date.get(Calendar.HOUR_OF_DAY); 200 int minute = date.get(Calendar.MINUTE); 201 int second = date.get(Calendar.SECOND); 202 203 double extra = (100.0 * year) + month - 190002.5; 204 double JD = (367.0 * year) 205 - (Math 206 .floor(7.0 * (year + Math.floor((month + 9.0) / 12.0)) / 4.0)) 207 + Math.floor((275.0 * month) / 9.0) + day 208 + ((hour + ((minute + (second / 60.0)) / 60.0)) / 24.0) 209 + 1721013.5 - ((0.5 * extra) / Math.abs(extra)) + 0.5; 210 return JD; 211 } 212 213 /** 214 * Convert Julian Day to centuries since J2000.0 215 * 216 * @param julian 217 * The Julian Day to convert 218 * @return the value corresponding to the Julian Day 219 */ 220 private static double julianDayToJulianCenturies(double julian) { 221 return (julian - 2451545) / 36525; 222 } 223 224 /** 225 * Convert centuries since J2000.0 to Julian Day 226 * 227 * @param t 228 * Number of Julian centuries since J2000.0 229 * @return The Julian Day corresponding to the value of t 230 */ 231 private static double julianCenturiesToJulianDay(double t) { 232 return (t * 36525) + 2451545; 233 } 234 235 /** 236 * Calculate the difference between true solar time and mean solar time 237 * 238 * @param t 239 * Number of Julian centuries since J2000.0 240 * @return 241 */ 242 private static double equationOfTime(double t) { 243 double epsilon = obliquityCorrection(t); 244 double l0 = geomMeanLongSun(t); 245 double e = eccentricityOfEarthsOrbit(t); 246 double m = geometricMeanAnomalyOfSun(t); 247 double y = Math.pow((Math.tan(Math.toRadians(epsilon) / 2)), 2); 248 249 double eTime = y * Math.sin(2 * Math.toRadians(l0)) - 2 * e 250 * Math.sin(Math.toRadians(m)) + 4 * e * y 251 * Math.sin(Math.toRadians(m)) 252 * Math.cos(2 * Math.toRadians(l0)) - 0.5 * y * y 253 * Math.sin(4 * Math.toRadians(l0)) - 1.25 * e * e 254 * Math.sin(2 * Math.toRadians(m)); 255 return Math.toDegrees(eTime) * 4; 256 } 257 258 /** 259 * Calculate the declination of the sun 260 * 261 * @param t 262 * Number of Julian centuries since J2000.0 263 * @return The Sun's declination in degrees 264 */ 265 private static double sunDeclination(double t) { 266 double e = obliquityCorrection(t); 267 double lambda = sunsApparentLongitude(t); 268 269 double sint = Math.sin(Math.toRadians(e)) 270 * Math.sin(Math.toRadians(lambda)); 271 return Math.toDegrees(Math.asin(sint)); 272 } 273 274 /** 275 * calculate the hour angle of the sun for a morning phenomenon for the 276 * given latitude 277 * 278 * @param lat 279 * Latitude of the observer in degrees 280 * @param solarDec 281 * declination of the sun in degrees 282 * @param zenithDistance 283 * zenith distance of the sun in degrees 284 * @return hour angle of sunrise in radians 285 */ 286 private static double hourAngleMorning(double lat, double solarDec, 287 double zenithDistance) { 288 return (Math.acos(Math.cos(Math.toRadians(zenithDistance)) 289 / (Math.cos(Math.toRadians(lat)) * Math.cos(Math 290 .toRadians(solarDec))) - Math.tan(Math.toRadians(lat)) 291 * Math.tan(Math.toRadians(solarDec)))); 292 } 293 294 /** 295 * Calculate the hour angle of the sun for an evening phenomenon for the 296 * given latitude 297 * 298 * @param lat 299 * Latitude of the observer in degrees 300 * @param solarDec 301 * declination of the Sun in degrees 302 * @param zenithDistance 303 * zenith distance of the sun in degrees 304 * @return hour angle of sunset in radians 305 */ 306 private static double hourAngleEvening(double lat, double solarDec, 307 double zenithDistance) { 308 return -hourAngleMorning(lat, solarDec, zenithDistance); 309 } 310 311 /** 312 * Calculate the corrected obliquity of the ecliptic 313 * 314 * @param t 315 * Number of Julian centuries since J2000.0 316 * @return Corrected obliquity in degrees 317 */ 318 private static double obliquityCorrection(double t) { 319 return meanObliquityOfEcliptic(t) + 0.00256 320 * Math.cos(Math.toRadians(125.04 - 1934.136 * t)); 321 } 322 323 /** 324 * Calculate the mean obliquity of the ecliptic 325 * 326 * @param t 327 * Number of Julian centuries since J2000.0 328 * @return Mean obliquity in degrees 329 */ 330 private static double meanObliquityOfEcliptic(double t) { 331 return 23 + (26 + (21.448 - t 332 * (46.815 + t * (0.00059 - t * (0.001813))) / 60)) / 60; 333 } 334 335 /** 336 * Calculate the geometric mean longitude of the sun 337 * 338 * @param t 339 * number of Julian centuries since J2000.0 340 * @return the geometric mean longitude of the sun in degrees 341 */ 342 private static double geomMeanLongSun(double t) { 343 double l0 = 280.46646 + t * (36000.76983 + 0.0003032 * t); 344 345 while ((l0 >= 0) && (l0 <= 360)) { 346 if (l0 > 360) { 347 l0 = l0 - 360; 348 } 349 350 if (l0 < 0) { 351 l0 = l0 + 360; 352 } 353 } 354 return l0; 355 } 356 357 /** 358 * Calculate the eccentricity of Earth's orbit 359 * 360 * @param t 361 * Number of Julian centuries since J2000.0 362 * @return the eccentricity 363 */ 364 private static double eccentricityOfEarthsOrbit(double t) { 365 return 0.016708634 - t * (0.000042037 + 0.0000001267 * t); 366 } 367 368 /** 369 * Calculate the geometric mean anomaly of the Sun 370 * 371 * @param t 372 * Number of Julian centuries since J2000.0 373 * @return the geometric mean anomaly of the Sun in degrees 374 */ 375 private static double geometricMeanAnomalyOfSun(double t) { 376 return 357.52911 + t * (35999.05029 - 0.0001537 * t); 377 } 378 379 /** 380 * Calculate the apparent longitude of the sun 381 * 382 * @param t 383 * Number of Julian centuries since J2000.0 384 * @return The apparent longitude of the Sun in degrees 385 */ 386 private static double sunsApparentLongitude(double t) { 387 return sunsTrueLongitude(t) - 0.00569 - 0.00478 388 * Math.sin(Math.toRadians(125.04 - 1934.136 * t)); 389 } 390 391 /** 392 * Calculate the true longitude of the sun 393 * 394 * @param t 395 * Number of Julian centuries since J2000.0 396 * @return The Sun's true longitude in degrees 397 */ 398 private static double sunsTrueLongitude(double t) { 399 return geomMeanLongSun(t) + equationOfCentreForSun(t); 400 } 401 402 /** 403 * Calculate the equation of centre for the Sun 404 * 405 * @param centuries 406 * Number of Julian centuries since J2000.0 407 * @return The equation of centre for the Sun in degrees 408 */ 409 private static double equationOfCentreForSun(double t) { 410 double m = geometricMeanAnomalyOfSun(t); 411 412 return Math.sin(Math.toRadians(m)) 413 * (1.914602 - t * (0.004817 + 0.000014 * t)) 414 + Math.sin(2 * Math.toRadians(m)) * (0.019993 - 0.000101 * t) 415 + Math.sin(3 * Math.toRadians(m)) * 0.000289; 416 } 417 }