001/* 002 * Zmanim Java API 003 * Copyright (C) 2004-2011 Eliyahu Hershfeld 004 * 005 * This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General 006 * Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) 007 * any later version. 008 * 009 * This library is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY; without even the implied 010 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more 011 * details. 012 * You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to 013 * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA, 014 * or connect to: http://www.gnu.org/licenses/old-licenses/lgpl-2.1.html 015 */ 016package net.sourceforge.zmanim.util; 017 018import java.util.Calendar; 019 020/** 021 * Implementation of sunrise and sunset methods to calculate astronomical times based on the <a 022 * href=""http://noaa.gov">NOAA</a> algorithm. This calculator uses the Java algorithm based on the implementation by <a 023 * href=""http://noaa.gov">NOAA - National Oceanic and Atmospheric Administration</a>'s <a href = 024 * "http://www.srrb.noaa.gov/highlights/sunrise/sunrise.html">Surface Radiation Research Branch</a>. NOAA's <a 025 * href="http://www.srrb.noaa.gov/highlights/sunrise/solareqns.PDF">implementation</a> is based on equations from <a 026 * href="http://www.willbell.com/math/mc1.htm">Astronomical Algorithms</a> by <a 027 * href="http://en.wikipedia.org/wiki/Jean_Meeus">Jean Meeus</a>. Added to the algorithm is an adjustment of the zenith 028 * to account for elevation. The algorithm can be found in the <a 029 * href="http://en.wikipedia.org/wiki/Sunrise_equation">Wikipedia Sunrise Equation</a> article. 030 * 031 * @author © Eliyahu Hershfeld 2011 032 * @version 0.1 033 */ 034public class NOAACalculator extends AstronomicalCalculator { 035 /** 036 * The <a href="http://en.wikipedia.org/wiki/Julian_day">Julian day</a> of January 1, 2000 037 */ 038 private static final double JULIAN_DAY_JAN_1_2000 = 2451545.0; 039 040 /** 041 * Julian days per century 042 */ 043 private static final double JULIAN_DAYS_PER_CENTURY = 36525.0; 044 045 /** 046 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getCalculatorName() 047 */ 048 public String getCalculatorName() { 049 return "US National Oceanic and Atmospheric Administration Algorithm"; 050 } 051 052 /** 053 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunrise(Calendar, GeoLocation, double, boolean) 054 */ 055 public double getUTCSunrise(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) { 056 double elevation = adjustForElevation ? geoLocation.getElevation() : 0; 057 double adjustedZenith = adjustZenith(zenith, elevation); 058 059 double sunrise = getSunriseUTC(getJulianDay(calendar), geoLocation.getLatitude(), -geoLocation.getLongitude(), 060 adjustedZenith); 061 sunrise = sunrise / 60; 062 063 // ensure that the time is >= 0 and < 24 064 while (sunrise < 0.0) { 065 sunrise += 24.0; 066 } 067 while (sunrise >= 24.0) { 068 sunrise -= 24.0; 069 } 070 return sunrise; 071 } 072 073 /** 074 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunset(Calendar, GeoLocation, double, boolean) 075 */ 076 public double getUTCSunset(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) { 077 double elevation = adjustForElevation ? geoLocation.getElevation() : 0; 078 double adjustedZenith = adjustZenith(zenith, elevation); 079 080 double sunset = getSunsetUTC(getJulianDay(calendar), geoLocation.getLatitude(), -geoLocation.getLongitude(), 081 adjustedZenith); 082 sunset = sunset / 60; 083 084 // ensure that the time is >= 0 and < 24 085 while (sunset < 0.0) { 086 sunset += 24.0; 087 } 088 while (sunset >= 24.0) { 089 sunset -= 24.0; 090 } 091 return sunset; 092 } 093 094 /** 095 * Return the <a href="http://en.wikipedia.org/wiki/Julian_day">Julian day</a> from a Java Calendar 096 * 097 * @param calendar 098 * The Java Calendar 099 * @return the Julian day corresponding to the date Note: Number is returned for start of day. Fractional days 100 * should be added later. 101 */ 102 private static double getJulianDay(Calendar calendar) { 103 int year = calendar.get(Calendar.YEAR); 104 int month = calendar.get(Calendar.MONTH) + 1; 105 int day = calendar.get(Calendar.DAY_OF_MONTH); 106 if (month <= 2) { 107 year -= 1; 108 month += 12; 109 } 110 int a = year / 100; 111 int b = 2 - a + a / 4; 112 113 return Math.floor(365.25 * (year + 4716)) + Math.floor(30.6001 * (month + 1)) + day + b - 1524.5; 114 } 115 116 /** 117 * Convert <a href="http://en.wikipedia.org/wiki/Julian_day">Julian day</a> to centuries since J2000.0. 118 * 119 * @param julianDay 120 * the Julian Day to convert 121 * @return the centuries since 2000 Julian corresponding to the Julian Day 122 */ 123 private static double getJulianCenturiesFromJulianDay(double julianDay) { 124 return (julianDay - JULIAN_DAY_JAN_1_2000) / JULIAN_DAYS_PER_CENTURY; 125 } 126 127 /** 128 * Convert centuries since J2000.0 to <a href="http://en.wikipedia.org/wiki/Julian_day">Julian day</a>. 129 * 130 * @param julianCenturies 131 * the number of Julian centuries since J2000.0 132 * @return the Julian Day corresponding to the Julian centuries passed in 133 */ 134 private static double getJulianDayFromJulianCenturies(double julianCenturies) { 135 return julianCenturies * JULIAN_DAYS_PER_CENTURY + JULIAN_DAY_JAN_1_2000; 136 } 137 138 /** 139 * Returns the Geometric <a href="http://en.wikipedia.org/wiki/Mean_longitude">Mean Longitude</a> of the Sun. 140 * 141 * @param julianCenturies 142 * the number of Julian centuries since J2000.0 143 * @return the Geometric Mean Longitude of the Sun in degrees 144 */ 145 private static double getSunGeometricMeanLongitude(double julianCenturies) { 146 double longitude = 280.46646 + julianCenturies * (36000.76983 + 0.0003032 * julianCenturies); 147 while (longitude > 360.0) { 148 longitude -= 360.0; 149 } 150 while (longitude < 0.0) { 151 longitude += 360.0; 152 } 153 154 return longitude; // in degrees 155 } 156 157 /** 158 * Returns the Geometric <a href="http://en.wikipedia.org/wiki/Mean_anomaly">Mean Anomaly</a> of the Sun. 159 * 160 * @param julianCenturies 161 * the number of Julian centuries since J2000.0 162 * @return the Geometric Mean Anomaly of the Sun in degrees 163 */ 164 private static double getSunGeometricMeanAnomaly(double julianCenturies) { 165 return 357.52911 + julianCenturies * (35999.05029 - 0.0001537 * julianCenturies); // in degrees 166 } 167 168 /** 169 * Return the <a href="http://en.wikipedia.org/wiki/Eccentricity_%28orbit%29">eccentricity of earth's orbit</a>. 170 * 171 * @param julianCenturies 172 * the number of Julian centuries since J2000.0 173 * @return the unitless eccentricity 174 */ 175 private static double getEarthOrbitEccentricity(double julianCenturies) { 176 return 0.016708634 - julianCenturies * (0.000042037 + 0.0000001267 * julianCenturies); // unitless 177 } 178 179 /** 180 * Returns the <a href="http://en.wikipedia.org/wiki/Equation_of_the_center">equation of center</a> for the sun. 181 * 182 * @param julianCenturies 183 * the number of Julian centuries since J2000.0 184 * @return the equation of center for the sun in degrees 185 */ 186 private static double getSunEquationOfCenter(double julianCenturies) { 187 double m = getSunGeometricMeanAnomaly(julianCenturies); 188 189 double mrad = Math.toRadians(m); 190 double sinm = Math.sin(mrad); 191 double sin2m = Math.sin(mrad + mrad); 192 double sin3m = Math.sin(mrad + mrad + mrad); 193 194 return sinm * (1.914602 - julianCenturies * (0.004817 + 0.000014 * julianCenturies)) + sin2m 195 * (0.019993 - 0.000101 * julianCenturies) + sin3m * 0.000289;// in degrees 196 } 197 198 /** 199 * Return the true longitude of the sun 200 * 201 * @param julianCenturies 202 * the number of Julian centuries since J2000.0 203 * @return the sun's true longitude in degrees 204 */ 205 private static double getSunTrueLongitude(double julianCenturies) { 206 double sunLongitude = getSunGeometricMeanLongitude(julianCenturies); 207 double center = getSunEquationOfCenter(julianCenturies); 208 209 return sunLongitude + center; // in degrees 210 } 211 212 // /** 213 // * Returns the <a href="http://en.wikipedia.org/wiki/True_anomaly">true anamoly</a> of the sun. 214 // * 215 // * @param julianCenturies 216 // * the number of Julian centuries since J2000.0 217 // * @return the sun's true anamoly in degrees 218 // */ 219 // private static double getSunTrueAnomaly(double julianCenturies) { 220 // double meanAnomaly = getSunGeometricMeanAnomaly(julianCenturies); 221 // double equationOfCenter = getSunEquationOfCenter(julianCenturies); 222 // 223 // return meanAnomaly + equationOfCenter; // in degrees 224 // } 225 226 /** 227 * Return the apparent longitude of the sun 228 * 229 * @param julianCenturies 230 * the number of Julian centuries since J2000.0 231 * @return sun's apparent longitude in degrees 232 */ 233 private static double getSunApparentLongitude(double julianCenturies) { 234 double sunTrueLongitude = getSunTrueLongitude(julianCenturies); 235 236 double omega = 125.04 - 1934.136 * julianCenturies; 237 double lambda = sunTrueLongitude - 0.00569 - 0.00478 * Math.sin(Math.toRadians(omega)); 238 return lambda; // in degrees 239 } 240 241 /** 242 * Returns the mean <a href="http://en.wikipedia.org/wiki/Axial_tilt">obliquity of the ecliptic</a> (Axial tilt). 243 * 244 * @param julianCenturies 245 * the number of Julian centuries since J2000.0 246 * @return the mean obliquity in degrees 247 */ 248 private static double getMeanObliquityOfEcliptic(double julianCenturies) { 249 double seconds = 21.448 - julianCenturies 250 * (46.8150 + julianCenturies * (0.00059 - julianCenturies * (0.001813))); 251 return 23.0 + (26.0 + (seconds / 60.0)) / 60.0; // in degrees 252 } 253 254 /** 255 * Returns the corrected <a href="http://en.wikipedia.org/wiki/Axial_tilt">obliquity of the ecliptic</a> (Axial 256 * tilt). 257 * 258 * @param julianCenturies 259 * the number of Julian centuries since J2000.0 260 * @return the corrected obliquity in degrees 261 */ 262 private static double getObliquityCorrection(double julianCenturies) { 263 double obliquityOfEcliptic = getMeanObliquityOfEcliptic(julianCenturies); 264 265 double omega = 125.04 - 1934.136 * julianCenturies; 266 return obliquityOfEcliptic + 0.00256 * Math.cos(Math.toRadians(omega)); // in degrees 267 } 268 269 /** 270 * Return the <a href="http://en.wikipedia.org/wiki/Declination">declination</a> of the sun. 271 * 272 * @param julianCenturies 273 * the number of Julian centuries since J2000.0 274 * @param sun 275 * 's declination in degrees 276 */ 277 private static double getSunDeclination(double julianCenturies) { 278 double obliquityCorrection = getObliquityCorrection(julianCenturies); 279 double lambda = getSunApparentLongitude(julianCenturies); 280 281 double sint = Math.sin(Math.toRadians(obliquityCorrection)) * Math.sin(Math.toRadians(lambda)); 282 double theta = Math.toDegrees(Math.asin(sint)); 283 return theta; // in degrees 284 } 285 286 /** 287 * Return the <a href="http://en.wikipedia.org/wiki/Equation_of_time">Equation of Time</a> - the difference between 288 * true solar time and mean solar time 289 * 290 * @param julianCenturies 291 * the number of Julian centuries since J2000.0 292 * @return equation of time in minutes of time 293 */ 294 private static double getEquationOfTime(double julianCenturies) { 295 double epsilon = getObliquityCorrection(julianCenturies); 296 double geomMeanLongSun = getSunGeometricMeanLongitude(julianCenturies); 297 double eccentricityEarthOrbit = getEarthOrbitEccentricity(julianCenturies); 298 double geomMeanAnomalySun = getSunGeometricMeanAnomaly(julianCenturies); 299 300 double y = Math.tan(Math.toRadians(epsilon) / 2.0); 301 y *= y; 302 303 double sin2l0 = Math.sin(2.0 * Math.toRadians(geomMeanLongSun)); 304 double sinm = Math.sin(Math.toRadians(geomMeanAnomalySun)); 305 double cos2l0 = Math.cos(2.0 * Math.toRadians(geomMeanLongSun)); 306 double sin4l0 = Math.sin(4.0 * Math.toRadians(geomMeanLongSun)); 307 double sin2m = Math.sin(2.0 * Math.toRadians(geomMeanAnomalySun)); 308 309 double equationOfTime = y * sin2l0 - 2.0 * eccentricityEarthOrbit * sinm + 4.0 * eccentricityEarthOrbit * y 310 * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * eccentricityEarthOrbit * eccentricityEarthOrbit * sin2m; 311 return Math.toDegrees(equationOfTime) * 4.0; // in minutes of time 312 } 313 314 /** 315 * Return the <a href="http://en.wikipedia.org/wiki/Hour_angle">hour angle</a> of the sun at sunrise for the 316 * latitude. 317 * 318 * @param lat 319 * , the latitude of observer in degrees 320 * @param solarDec 321 * the declination angle of sun in degrees 322 * @return hour angle of sunrise in radians 323 */ 324 private static double getSunHourAngleAtSunrise(double lat, double solarDec, double zenith) { 325 double latRad = Math.toRadians(lat); 326 double sdRad = Math.toRadians(solarDec); 327 328 return (Math.acos(Math.cos(Math.toRadians(zenith)) / (Math.cos(latRad) * Math.cos(sdRad)) - Math.tan(latRad) 329 * Math.tan(sdRad))); // in radians 330 } 331 332 /** 333 * Returns the <a href="http://en.wikipedia.org/wiki/Hour_angle">hour angle</a> of the sun at sunset for the 334 * latitude. TODO: use - {@link #getSunHourAngleAtSunrise(double, double, double)} implementation to avoid 335 * duplication of code. 336 * 337 * @param lat 338 * the latitude of observer in degrees 339 * @param solarDec 340 * the declination angle of sun in degrees 341 * @return the hour angle of sunset in radians 342 */ 343 private static double getSunHourAngleAtSunset(double lat, double solarDec, double zenith) { 344 double latRad = Math.toRadians(lat); 345 double sdRad = Math.toRadians(solarDec); 346 347 double hourAngle = (Math.acos(Math.cos(Math.toRadians(zenith)) / (Math.cos(latRad) * Math.cos(sdRad)) 348 - Math.tan(latRad) * Math.tan(sdRad))); 349 return -hourAngle; // in radians 350 } 351 352 /** 353 * Return the <a href="http://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC) 354 * of sunrise for the given day at the given location on earth 355 * 356 * @param julianDay 357 * the Julian day 358 * @param latitude 359 * the latitude of observer in degrees 360 * @param longitude 361 * the longitude of observer in degrees 362 * @return the time in minutes from zero UTC 363 */ 364 private static double getSunriseUTC(double julianDay, double latitude, double longitude, double zenith) { 365 double julianCenturies = getJulianCenturiesFromJulianDay(julianDay); 366 367 // Find the time of solar noon at the location, and use that declination. This is better than start of the 368 // Julian day 369 370 double noonmin = getSolarNoonUTC(julianCenturies, longitude); 371 double tnoon = getJulianCenturiesFromJulianDay(julianDay + noonmin / 1440.0); 372 373 // First pass to approximate sunrise (using solar noon) 374 375 double eqTime = getEquationOfTime(tnoon); 376 double solarDec = getSunDeclination(tnoon); 377 double hourAngle = getSunHourAngleAtSunrise(latitude, solarDec, zenith); 378 379 double delta = longitude - Math.toDegrees(hourAngle); 380 double timeDiff = 4 * delta; // in minutes of time 381 double timeUTC = 720 + timeDiff - eqTime; // in minutes 382 383 // Second pass includes fractional Julian Day in gamma calc 384 385 double newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + timeUTC 386 / 1440.0); 387 eqTime = getEquationOfTime(newt); 388 solarDec = getSunDeclination(newt); 389 hourAngle = getSunHourAngleAtSunrise(latitude, solarDec, zenith); 390 delta = longitude - Math.toDegrees(hourAngle); 391 timeDiff = 4 * delta; 392 timeUTC = 720 + timeDiff - eqTime; // in minutes 393 return timeUTC; 394 } 395 396 /** 397 * Return the <a href="http://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC) 398 * of <a href="http://en.wikipedia.org/wiki/Noon#Solar_noon">solar noon</a> for the given day at the given location 399 * on earth. 400 * 401 * @param julianCenturies 402 * the number of Julian centuries since J2000.0 403 * @param longitude 404 * the longitude of observer in degrees 405 * @return the time in minutes from zero UTC 406 */ 407 private static double getSolarNoonUTC(double julianCenturies, double longitude) { 408 // First pass uses approximate solar noon to calculate eqtime 409 double tnoon = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + longitude 410 / 360.0); 411 double eqTime = getEquationOfTime(tnoon); 412 double solNoonUTC = 720 + (longitude * 4) - eqTime; // min 413 414 double newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) - 0.5 415 + solNoonUTC / 1440.0); 416 417 eqTime = getEquationOfTime(newt); 418 return 720 + (longitude * 4) - eqTime; // min 419 } 420 421 /** 422 * Return the <a href="http://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC) 423 * of sunset for the given day at the given location on earth 424 * 425 * @param julianDay 426 * the Julian day 427 * @param latitude 428 * the latitude of observer in degrees 429 * @param longitude 430 * : longitude of observer in degrees 431 * @param zenith 432 * @return the time in minutes from zero Universal Coordinated Time (UTC) 433 */ 434 private static double getSunsetUTC(double julianDay, double latitude, double longitude, double zenith) { 435 double julianCenturies = getJulianCenturiesFromJulianDay(julianDay); 436 437 // Find the time of solar noon at the location, and use that declination. This is better than start of the 438 // Julian day 439 440 double noonmin = getSolarNoonUTC(julianCenturies, longitude); 441 double tnoon = getJulianCenturiesFromJulianDay(julianDay + noonmin / 1440.0); 442 443 // First calculates sunrise and approx length of day 444 445 double eqTime = getEquationOfTime(tnoon); 446 double solarDec = getSunDeclination(tnoon); 447 double hourAngle = getSunHourAngleAtSunset(latitude, solarDec, zenith); 448 449 double delta = longitude - Math.toDegrees(hourAngle); 450 double timeDiff = 4 * delta; 451 double timeUTC = 720 + timeDiff - eqTime; 452 453 // Second pass includes fractional Julian Day in gamma calc 454 455 double newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + timeUTC 456 / 1440.0); 457 eqTime = getEquationOfTime(newt); 458 solarDec = getSunDeclination(newt); 459 hourAngle = getSunHourAngleAtSunset(latitude, solarDec, zenith); 460 461 delta = longitude - Math.toDegrees(hourAngle); 462 timeDiff = 4 * delta; 463 timeUTC = 720 + timeDiff - eqTime; // in minutes 464 return timeUTC; 465 } 466}