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 }