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