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 }