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