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 }