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 based on the <a href=""http://noaa.gov">NOAA</a> algorithm. 024 * This calculator uses the Java algorithm based on the implementation by <a 025 * href=""http://noaa.gov">NOAA - National Oceanic and Atmospheric 026 * Administration</a>'s <a href = 027 * "http://www.srrb.noaa.gov/highlights/sunrise/sunrisehtml">Surface Radiation 028 * Research Branch</a>. NOAA's <a 029 * href="http://www.srrb.noaa.gov/highlights/sunrise/solareqns.PDF">implementation</a> 030 * is based on equations from <a 031 * href="http://www.willbell.com/math/mc1.htm">Astronomical Algorithms</a> by 032 * <a href="http://en.wikipedia.org/wiki/Jean_Meeus">Jean Meeus</a>. Added to 033 * the algorithm is an adjustment of the zenith to account for elevation. 034 * 035 * @author © Eliyahu Hershfeld 2011 036 * @version 0.1 037 */ 038 public class NOAACalculator extends AstronomicalCalculator { 039 private String calculatorName = "US National Oceanic and Atmospheric Administration Algorithm"; 040 public String getCalculatorName() { 041 return this.calculatorName; 042 } 043 044 /** 045 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunrise(AstronomicalCalendar, 046 * double, boolean) 047 */ 048 public double getUTCSunrise(AstronomicalCalendar astronomicalCalendar, 049 double zenith, boolean adjustForElevation) { 050 double adjustedZenith = zenith; 051 052 // if (astronomicalCalendar.getCalendar().get(Calendar.YEAR) <= 2000) { 053 // throw new ZmanimException( 054 // "NOAACalculator can not calculate times earlier than the year 2000. Please try a date with a different year."); 055 // } 056 057 if (adjustForElevation) { 058 adjustedZenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation() 059 .getElevation()); 060 } else { 061 adjustedZenith = adjustZenith(zenith, 0); 062 } 063 064 double sunRise = calcSunriseUTC(calcJD(astronomicalCalendar 065 .getCalendar()), astronomicalCalendar.getGeoLocation() 066 .getLatitude(), -astronomicalCalendar.getGeoLocation() 067 .getLongitude(), adjustedZenith); 068 return sunRise / 60; 069 } 070 071 /** 072 * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunset(AstronomicalCalendar, 073 * double, boolean) 074 * @throws ZmanimException 075 * if the year entered < 2001. This calculator can't properly 076 * deal with the year 2000. It can properly calculate times for 077 * years <> 2000. 078 */ 079 public double getUTCSunset(AstronomicalCalendar astronomicalCalendar, 080 double zenith, boolean adjustForElevation) { 081 double adjustedZenith = zenith; 082 if (adjustForElevation) { 083 adjustedZenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation() 084 .getElevation()); 085 } else { 086 adjustedZenith = adjustZenith(zenith, 0); 087 } 088 089 double sunSet = calcSunsetUTC( 090 calcJD(astronomicalCalendar.getCalendar()), 091 astronomicalCalendar.getGeoLocation().getLatitude(), 092 -astronomicalCalendar.getGeoLocation().getLongitude(), adjustedZenith); 093 return sunSet / 60; 094 } 095 096 /** 097 * Generate a Julian day from Java Calendar 098 * 099 * @param date 100 * Java Calendar 101 * @return the Julian day corresponding to the date Note: Number is returned 102 * for start of day. Fractional days should be added later. 103 */ 104 private static double calcJD(Calendar date) { 105 int year = date.get(Calendar.YEAR); 106 int month = date.get(Calendar.MONTH) + 1; 107 int day = date.get(Calendar.DAY_OF_MONTH); 108 if (month <= 2) { 109 year -= 1; 110 month += 12; 111 } 112 double A = Math.floor(year / 100); 113 double B = 2 - A + Math.floor(A / 4); 114 115 return Math.floor(365.25 * (year + 4716)) 116 + Math.floor(30.6001 * (month + 1)) + day + B - 1524.5; 117 } 118 119 /** 120 * convert Julian Day to centuries since J2000.0. 121 * 122 * @param jd 123 * the Julian Day to convert 124 * @return the T value corresponding to the Julian Day 125 */ 126 private static double calcTimeJulianCent(double jd) { 127 return (jd - 2451545.0) / 36525.0; 128 } 129 130 /** 131 * Convert centuries since J2000.0 to Julian Day. 132 * 133 * @param t 134 * the number of Julian centuries since J2000.0 135 * @return the Julian Day corresponding to the t value 136 */ 137 private static double calcJDFromJulianCent(double t) { 138 return t * 36525.0 + 2451545.0; 139 } 140 141 /** 142 * calculates the Geometric Mean Longitude of the Sun 143 * 144 * @param t 145 * the number of Julian centuries since J2000.0 146 * @return the Geometric Mean Longitude of the Sun in degrees 147 */ 148 private static double calcGeomMeanLongSun(double t) { 149 double L0 = 280.46646 + t * (36000.76983 + 0.0003032 * t); 150 while (L0 > 360.0) { 151 L0 -= 360.0; 152 } 153 while (L0 < 0.0) { 154 L0 += 360.0; 155 } 156 157 return L0; // in degrees 158 } 159 160 /** 161 * Calculate the Geometric Mean Anomaly of the Sun 162 * 163 * @param t 164 * the number of Julian centuries since J2000.0 165 * @return the Geometric Mean Anomaly of the Sun in degrees 166 */ 167 private static double calcGeomMeanAnomalySun(double t) { 168 double M = 357.52911 + t * (35999.05029 - 0.0001537 * t); 169 return M; // in degrees 170 } 171 172 /** 173 * calculate the eccentricity of earth's orbit 174 * 175 * @param t 176 * the number of Julian centuries since J2000.0 177 * @return the unitless eccentricity 178 */ 179 private static double calcEccentricityEarthOrbit(double t) { 180 double e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t); 181 return e; // unitless 182 } 183 184 /** 185 * Calculate the equation of center for the sun 186 * 187 * @param t 188 * the number of Julian centuries since J2000.0 189 * @return the equation of center for the sun in degrees 190 */ 191 private static double calcSunEqOfCenter(double t) { 192 double m = calcGeomMeanAnomalySun(t); 193 194 double mrad = Math.toRadians(m); 195 double sinm = Math.sin(mrad); 196 double sin2m = Math.sin(mrad + mrad); 197 double sin3m = Math.sin(mrad + mrad + mrad); 198 199 double C = sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) + sin2m 200 * (0.019993 - 0.000101 * t) + sin3m * 0.000289; 201 return C; // in degrees 202 } 203 204 /** 205 * Calculate the true longitude of the sun 206 * 207 * @param t 208 * the number of Julian centuries since J2000.0 209 * @return the sun's true longitude in degrees 210 */ 211 private static double calcSunTrueLong(double t) { 212 double l0 = calcGeomMeanLongSun(t); 213 double c = calcSunEqOfCenter(t); 214 215 double O = l0 + c; 216 return O; // in degrees 217 } 218 219 // /** 220 // * Calculate the true anamoly of the sun 221 // * 222 // * @param t 223 // * the number of Julian centuries since J2000.0 224 // * @return the sun's true anamoly in degrees 225 // */ 226 // private static double calcSunTrueAnomaly(double t) { 227 // double m = calcGeomMeanAnomalySun(t); 228 // double c = calcSunEqOfCenter(t); 229 // 230 // double v = m + c; 231 // return v; // in degrees 232 // } 233 234 /** 235 * calculate the apparent longitude of the sun 236 * 237 * @param t 238 * the number of Julian centuries since J2000.0 239 * @return sun's apparent longitude in degrees 240 */ 241 private static double calcSunApparentLong(double t) { 242 double o = calcSunTrueLong(t); 243 244 double omega = 125.04 - 1934.136 * t; 245 double lambda = o - 0.00569 - 0.00478 * Math.sin(Math.toRadians(omega)); 246 return lambda; // in degrees 247 } 248 249 /** 250 * Calculate the mean obliquity of the ecliptic 251 * 252 * @param t 253 * the number of Julian centuries since J2000.0 254 * @return the mean obliquity in degrees 255 */ 256 private static double calcMeanObliquityOfEcliptic(double t) { 257 double seconds = 21.448 - t 258 * (46.8150 + t * (0.00059 - t * (0.001813))); 259 double e0 = 23.0 + (26.0 + (seconds / 60.0)) / 60.0; 260 return e0; // in degrees 261 } 262 263 /** 264 * calculate the corrected obliquity of the ecliptic 265 * 266 * @param t 267 * the number of Julian centuries since J2000.0 268 * @return the corrected obliquity in degrees 269 */ 270 private static double calcObliquityCorrection(double t) { 271 double e0 = calcMeanObliquityOfEcliptic(t); 272 273 double omega = 125.04 - 1934.136 * t; 274 double e = e0 + 0.00256 * Math.cos(Math.toRadians(omega)); 275 return e; // in degrees 276 } 277 278 /** 279 * Calculate the declination of the sun 280 * 281 * @param t 282 * the number of Julian centuries since J2000.0 283 * @param sun's 284 * declination in degrees 285 */ 286 private static double calcSunDeclination(double t) { 287 double e = calcObliquityCorrection(t); 288 double lambda = calcSunApparentLong(t); 289 290 double sint = Math.sin(Math.toRadians(e)) 291 * Math.sin(Math.toRadians(lambda)); 292 double theta = Math.toDegrees(Math.asin(sint)); 293 return theta; // in degrees 294 } 295 296 /** 297 * calculate the difference between true solar time and mean solar time 298 * 299 * @param t 300 * the number of Julian centuries since J2000.0 301 * @return equation of time in minutes of time 302 */ 303 private static double calcEquationOfTime(double t) { 304 double epsilon = calcObliquityCorrection(t); 305 double l0 = calcGeomMeanLongSun(t); 306 double e = calcEccentricityEarthOrbit(t); 307 double m = calcGeomMeanAnomalySun(t); 308 309 double y = Math.tan(Math.toRadians(epsilon) / 2.0); 310 y *= y; 311 312 double sin2l0 = Math.sin(2.0 * Math.toRadians(l0)); 313 double sinm = Math.sin(Math.toRadians(m)); 314 double cos2l0 = Math.cos(2.0 * Math.toRadians(l0)); 315 double sin4l0 = Math.sin(4.0 * Math.toRadians(l0)); 316 double sin2m = Math.sin(2.0 * Math.toRadians(m)); 317 318 double Etime = y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm 319 * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m; 320 return Math.toDegrees(Etime) * 4.0; // in minutes of time 321 } 322 323 /** 324 * Calculate the hour angle of the sun at sunrise for the latitude 325 * 326 * @param lat, 327 * the latitude of observer in degrees 328 * @param solarDec 329 * the declination angle of sun in degrees 330 * @return hour angle of sunrise in radians 331 */ 332 private static double calcHourAngleSunrise(double lat, double solarDec, 333 double zenith) { 334 double latRad = Math.toRadians(lat); 335 double sdRad = Math.toRadians(solarDec); 336 337 // double HAarg = 338 // (Math.cos(Math.toRadians(zenith))/(Math.cos(latRad)*Math.cos(sdRad))-Math.tan(latRad) 339 // * Math.tan(sdRad)); 340 341 double HA = (Math.acos(Math.cos(Math.toRadians(zenith)) 342 / (Math.cos(latRad) * Math.cos(sdRad)) - Math.tan(latRad) 343 * Math.tan(sdRad))); 344 return HA; // in radians 345 } 346 347 /** 348 * Calculate the hour angle of the sun at sunset for the latitude 349 * 350 * @param lat 351 * the latitude of observer in degrees 352 * @param solarDec 353 * the declination angle of sun in degrees 354 * @return the hour angle of sunset in radians 355 * TODO: use - calcHourAngleSunrise implementation 356 */ 357 private static double calcHourAngleSunset(double lat, double solarDec, 358 double zenith) { 359 double latRad = Math.toRadians(lat); 360 double sdRad = Math.toRadians(solarDec); 361 362 // double HAarg = 363 // (Math.cos(Math.toRadians(zenith))/(Math.cos(latRad)*Math.cos(sdRad))-Math.tan(latRad) 364 // * Math.tan(sdRad)); 365 366 double HA = (Math.acos(Math.cos(Math.toRadians(zenith)) 367 / (Math.cos(latRad) * Math.cos(sdRad)) - Math.tan(latRad) 368 * Math.tan(sdRad))); 369 return -HA; // in radians 370 } 371 372 /** 373 * Calculate the Universal Coordinated Time (UTC) of sunrise for the given 374 * day at the given location on earth 375 * 376 * @param JD 377 * the julian day 378 * @param latitude 379 * the latitude of observer in degrees 380 * @param longitude 381 * the longitude of observer in degrees 382 * @return the time in minutes from zero Z 383 */ 384 private static double calcSunriseUTC(double JD, double latitude, 385 double longitude, double zenith) { 386 double t = calcTimeJulianCent(JD); 387 388 // *** Find the time of solar noon at the location, and use 389 // that declination. This is better than start of the 390 // Julian day 391 392 double noonmin = calcSolNoonUTC(t, longitude); 393 double tnoon = calcTimeJulianCent(JD + noonmin / 1440.0); 394 395 // *** First pass to approximate sunrise (using solar noon) 396 397 double eqTime = calcEquationOfTime(tnoon); 398 double solarDec = calcSunDeclination(tnoon); 399 double hourAngle = calcHourAngleSunrise(latitude, solarDec, zenith); 400 401 double delta = longitude - Math.toDegrees(hourAngle); 402 double timeDiff = 4 * delta; // in minutes of time 403 double timeUTC = 720 + timeDiff - eqTime; // in minutes 404 405 // *** Second pass includes fractional jday in gamma calc 406 407 double newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC 408 / 1440.0); 409 eqTime = calcEquationOfTime(newt); 410 solarDec = calcSunDeclination(newt); 411 hourAngle = calcHourAngleSunrise(latitude, solarDec, zenith); 412 delta = longitude - Math.toDegrees(hourAngle); 413 timeDiff = 4 * delta; 414 timeUTC = 720 + timeDiff - eqTime; // in minutes 415 return timeUTC; 416 } 417 418 /** 419 * calculate the Universal Coordinated Time (UTC) of solar noon for the 420 * given day at the given location on earth 421 * 422 * @param t 423 * the number of Julian centuries since J2000.0 424 * @param longitude 425 * the longitude of observer in degrees 426 * @return the time in minutes from zero Z 427 */ 428 private static double calcSolNoonUTC(double t, double longitude) { 429 // First pass uses approximate solar noon to calculate eqtime 430 double tnoon = calcTimeJulianCent(calcJDFromJulianCent(t) + longitude 431 / 360.0); 432 double eqTime = calcEquationOfTime(tnoon); 433 double solNoonUTC = 720 + (longitude * 4) - eqTime; // min 434 435 double newt = calcTimeJulianCent(calcJDFromJulianCent(t) - 0.5 436 + solNoonUTC / 1440.0); 437 438 eqTime = calcEquationOfTime(newt); 439 return 720 + (longitude * 4) - eqTime; // min 440 } 441 442 /** 443 * calculate the Universal Coordinated Time (UTC) of sunset for the given 444 * day at the given location on earth 445 * 446 * @param JD 447 * the julian day 448 * @param latitude 449 * the latitude of observer in degrees 450 * @param longitude : 451 * longitude of observer in degrees 452 * @param zenith 453 * @return the time in minutes from zero Z 454 */ 455 private static double calcSunsetUTC(double JD, double latitude, 456 double longitude, double zenith) { 457 double t = calcTimeJulianCent(JD); 458 459 // *** Find the time of solar noon at the location, and use 460 // that declination. This is better than start of the 461 // Julian day 462 463 double noonmin = calcSolNoonUTC(t, longitude); 464 double tnoon = calcTimeJulianCent(JD + noonmin / 1440.0); 465 466 // First calculates sunrise and approx length of day 467 468 double eqTime = calcEquationOfTime(tnoon); 469 double solarDec = calcSunDeclination(tnoon); 470 double hourAngle = calcHourAngleSunset(latitude, solarDec, zenith); 471 472 double delta = longitude - Math.toDegrees(hourAngle); 473 double timeDiff = 4 * delta; 474 double timeUTC = 720 + timeDiff - eqTime; 475 476 // first pass used to include fractional day in gamma calc 477 478 double newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC 479 / 1440.0); 480 eqTime = calcEquationOfTime(newt); 481 solarDec = calcSunDeclination(newt); 482 hourAngle = calcHourAngleSunset(latitude, solarDec, zenith); 483 484 delta = longitude - Math.toDegrees(hourAngle); 485 timeDiff = 4 * delta; 486 return 720 + timeDiff - eqTime; // in minutes 487 } 488 }