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