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