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