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 &copy; 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    }