001 /* 002 * Zmanim Java API 003 * Copyright (C) 2004-2007 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 net.sourceforge.zmanim.ZmanimException; 021 import java.util.Calendar; 022 023 /** 024 * Implementation of sunrise and sunset methods to calculate astronomical times. 025 * This calculator uses the Java algorithm written by <a 026 * href="http://www.jstot.me.uk/jsuntimes/">Jonathan Stott</a> that is based on 027 * the implementation by <a href=""http://noaa.gov">NOAA - National Oceanic and 028 * Atmospheric Administration</a>'s <a href = 029 * "http://www.srrb.noaa.gov/highlights/sunrise/sunrisehtml">Surface Radiation 030 * Research Branch</a>. NOAA's <a 031 * href="http://www.srrb.noaa.gov/highlights/sunrise/solareqns.PDF">implementation</a> 032 * is based on equations from <a 033 * href="http://www.willbell.com/math/mc1.htm">Astronomical Algorithms</a> by 034 * <a href="http://en.wikipedia.org/wiki/Jean_Meeus">Jean Meeus</a>. Jonathan's 035 * implementation was released under the GPL. Added to the algorithm is an 036 * adjustment of the zenith to account for elevation. 037 * 038 * @author Jonathan Stott 2000 - 2004 039 * @author © Eliyahu Hershfeld 2004 - 2007 040 * @version 1.1 041 * @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. 042 * @see net.sourceforge.zmanim.util.NOAACalculator 043 */ 044 public class JSuntimeCalculator extends AstronomicalCalculator { 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 "US National Oceanic and Atmospheric Administration Algorithm"; 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 //if(timeMins/60 < 4) {System.out.println("zenith="+zenith);} 080 return timeMins / 60; 081 } 082 083 /** 084 * @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. 085 * @see net.sourceforge.zmanim.util.NOAACalculator#getUTCSunset(AstronomicalCalendar, double, boolean) 086 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunset(AstronomicalCalendar, 087 * double, boolean) 088 * @throws ZmanimException 089 * if the year entered == 2000. This calculator can't properly 090 * deal with the year 2000. It can properly calculate times for 091 * years <> 2000. 092 */ 093 public double getUTCSunset(AstronomicalCalendar astronomicalCalendar, double zenith, boolean adjustForElevation) { 094 if (astronomicalCalendar.getCalendar().get(Calendar.YEAR) == 2000) { 095 throw new ZmanimException( 096 "JSuntimeCalculator can not calculate times for the year 2000. Please try a date with a different year."); 097 } 098 099 if (adjustForElevation) { 100 zenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation() 101 .getElevation()); 102 } else { 103 zenith = adjustZenith(zenith, 0); 104 } 105 double timeMins = eveningPhenomenon(dateToJulian(astronomicalCalendar 106 .getCalendar()), astronomicalCalendar.getGeoLocation() 107 .getLatitude(), -astronomicalCalendar.getGeoLocation() 108 .getLongitude(), zenith); 109 //if(timeMins/60 > 28) {System.out.println("zenith="+zenith);} 110 return timeMins / 60; 111 } 112 113 /** 114 * Calculate the UTC of a morning phenomenon for the given day at the given 115 * latitude and longitude on Earth 116 * 117 * @param julian 118 * Julian day 119 * @param latitude 120 * latitude of observer in degrees 121 * @param longitude 122 * longitude of observer in degrees 123 * @param zenithDistance 124 * one of Sun.SUNRISE_SUNSET_ZENITH_DISTANCE, 125 * Sun.CIVIL_TWILIGHT_ZENITH_DISTANCE, 126 * Sun.NAUTICAL_TWILIGHT_ZENITH_DISTANCE, 127 * Sun.ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE. 128 * @return time in minutes from zero Z 129 */ 130 private static double morningPhenomenon(double julian, double latitude, 131 double longitude, double zenithDistance) { 132 double t = julianDayToJulianCenturies(julian); 133 double eqtime = equationOfTime(t); 134 double solarDec = sunDeclination(t); 135 double hourangle = hourAngleMorning(latitude, solarDec, zenithDistance); 136 double delta = longitude - Math.toDegrees(hourangle); 137 double timeDiff = 4 * delta; 138 double timeUTC = 720 + timeDiff - eqtime; 139 140 // Second pass includes fractional julian day in gamma calc 141 double newt = julianDayToJulianCenturies(julianCenturiesToJulianDay(t) 142 + timeUTC / 1440); 143 eqtime = equationOfTime(newt); 144 solarDec = sunDeclination(newt); 145 hourangle = hourAngleMorning(latitude, solarDec, zenithDistance); 146 delta = longitude - Math.toDegrees(hourangle); 147 timeDiff = 4 * delta; 148 149 return 720 + timeDiff - eqtime; 150 } 151 152 /** 153 * Calculate the UTC of an evening phenomenon for the given day at the given 154 * latitude and longitude on Earth 155 * 156 * @param julian 157 * Julian day 158 * @param latitude 159 * latitude of observer in degrees 160 * @param longitude 161 * longitude of observer in degrees 162 * @param zenithDistance 163 * one of Sun.SUNRISE_SUNSET_ZENITH_DISTANCE, 164 * Sun.CIVIL_TWILIGHT_ZENITH_DISTANCE, 165 * Sun.NAUTICAL_TWILIGHT_ZENITH_DISTANCE, 166 * Sun.ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE. 167 * @return time in minutes from zero Z 168 */ 169 private static double eveningPhenomenon(double julian, double latitude, 170 double longitude, double zenithDistance) { 171 double t = julianDayToJulianCenturies(julian); 172 173 // First calculates sunrise and approx length of day 174 double eqtime = equationOfTime(t); 175 double solarDec = sunDeclination(t); 176 double hourangle = hourAngleEvening(latitude, solarDec, zenithDistance); 177 178 double delta = longitude - Math.toDegrees(hourangle); 179 double timeDiff = 4 * delta; 180 double timeUTC = 720 + timeDiff - eqtime; 181 182 // first pass used to include fractional day in gamma calc 183 double newt = julianDayToJulianCenturies(julianCenturiesToJulianDay(t) 184 + timeUTC / 1440); 185 eqtime = equationOfTime(newt); 186 solarDec = sunDeclination(newt); 187 hourangle = hourAngleEvening(latitude, solarDec, zenithDistance); 188 189 delta = longitude - Math.toDegrees(hourangle); 190 timeDiff = 4 * delta; 191 192 return 720 + timeDiff - eqtime; 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 return (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 } 211 212 /** 213 * Convert Julian Day to centuries since J2000.0 214 * 215 * @param julian 216 * The Julian Day to convert 217 * @return the value corresponding to the Julian Day 218 */ 219 private static double julianDayToJulianCenturies(double julian) { 220 return (julian - 2451545) / 36525; 221 } 222 223 /** 224 * Convert centuries since J2000.0 to Julian Day 225 * 226 * @param t 227 * Number of Julian centuries since J2000.0 228 * @return The Julian Day corresponding to the value of t 229 */ 230 private static double julianCenturiesToJulianDay(double t) { 231 return (t * 36525) + 2451545; 232 } 233 234 /** 235 * Calculate the difference between true solar time and mean solar time 236 * 237 * @param t 238 * Number of Julian centuries since J2000.0 239 * @return 240 */ 241 private static double equationOfTime(double t) { 242 double epsilon = obliquityCorrection(t); 243 double l0 = geomMeanLongSun(t); 244 double e = eccentricityOfEarthsOrbit(t); 245 double m = geometricMeanAnomalyOfSun(t); 246 double y = Math.pow((Math.tan(Math.toRadians(epsilon) / 2)), 2); 247 248 double eTime = y * Math.sin(2 * Math.toRadians(l0)) - 2 * e 249 * Math.sin(Math.toRadians(m)) + 4 * e * y 250 * Math.sin(Math.toRadians(m)) 251 * Math.cos(2 * Math.toRadians(l0)) - 0.5 * y * y 252 * Math.sin(4 * Math.toRadians(l0)) - 1.25 * e * e 253 * Math.sin(2 * Math.toRadians(m)); 254 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 355 return l0; 356 } 357 358 /** 359 * Calculate the eccentricity of Earth's orbit 360 * 361 * @param t 362 * Number of Julian centuries since J2000.0 363 * @return the eccentricity 364 */ 365 private static double eccentricityOfEarthsOrbit(double t) { 366 return 0.016708634 - t * (0.000042037 + 0.0000001267 * t); 367 } 368 369 /** 370 * Calculate the geometric mean anomaly of the Sun 371 * 372 * @param t 373 * Number of Julian centuries since J2000.0 374 * @return the geometric mean anomaly of the Sun in degrees 375 */ 376 private static double geometricMeanAnomalyOfSun(double t) { 377 return 357.52911 + t * (35999.05029 - 0.0001537 * t); 378 } 379 380 /** 381 * Calculate the apparent longitude of the sun 382 * 383 * @param t 384 * Number of Julian centuries since J2000.0 385 * @return The apparent longitude of the Sun in degrees 386 */ 387 private static double sunsApparentLongitude(double t) { 388 return sunsTrueLongitude(t) - 0.00569 - 0.00478 389 * Math.sin(Math.toRadians(125.04 - 1934.136 * t)); 390 } 391 392 /** 393 * Calculate the true longitude of the sun 394 * 395 * @param t 396 * Number of Julian centuries since J2000.0 397 * @return The Sun's true longitude in degrees 398 */ 399 private static double sunsTrueLongitude(double t) { 400 return geomMeanLongSun(t) + equationOfCentreForSun(t); 401 } 402 403 /** 404 * Calculate the equation of centre for the Sun 405 * 406 * @param centuries 407 * Number of Julian centuries since J2000.0 408 * @return The equation of centre for the Sun in degrees 409 */ 410 private static double equationOfCentreForSun(double t) { 411 double m = geometricMeanAnomalyOfSun(t); 412 413 return Math.sin(Math.toRadians(m)) 414 * (1.914602 - t * (0.004817 + 0.000014 * t)) 415 + Math.sin(2 * Math.toRadians(m)) * (0.019993 - 0.000101 * t) 416 + Math.sin(3 * Math.toRadians(m)) * 0.000289; 417 } 418 }