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    }