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.
024     * This calculator uses the Java algorithm written by <a
025     * href="http://www.jstot.me.uk/jsuntimes/">Jonathan Stott</a> that is based on
026     * the implementation by <a href=""http://noaa.gov">NOAA - National Oceanic and
027     * Atmospheric 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>. Jonathan's
034     * implementation was released under the GPL. Added to the algorithm is an
035     * adjustment of the zenith to account for elevation.
036     *
037     * @author Jonathan Stott 2000 - 2004
038     * @author &copy; Eliyahu Hershfeld 2004 - 2008
039     * @version 1.1
040     * @deprecated This class is based on the NOAA algorithm but does not return calculations that match the NOAA algorithm JavaScript implementation. The calculations are about 2 minutes off. This call has been replaced by the NOAACalculator class.
041     * @see net.sourceforge.zmanim.util.NOAACalculator
042     */
043    public class JSuntimeCalculator extends AstronomicalCalculator {
044            private String calculatorName =  "US National Oceanic and Atmospheric Administration Algorithm";
045            /**
046             * @deprecated This class is based on the NOAA algorithm but does not return calculations that match the NOAA algorithm JavaScript implementation. The calculations are about 2 minutes off. This call has been replaced by the NOAACalculator class.
047             * @see net.sourceforge.zmanim.util.NOAACalculator#getCalculatorName()
048             */
049            public String getCalculatorName() {
050                    return calculatorName;
051            }
052    
053            /**
054             * @deprecated This class is based on the NOAA algorithm but does not return calculations that match the NOAA algorithm JavaScript implementation. The calculations are about 2 minutes off. This call has been replaced by the NOAACalculator class.
055             * @see net.sourceforge.zmanim.util.NOAACalculator#getUTCSunrise(AstronomicalCalendar, double, boolean)
056             * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunrise(AstronomicalCalendar,
057             *      double, boolean)
058             * @throws ZmanimException
059             *             if the year entered == 2000. This calculator can't properly
060             *             deal with the year 2000. It can properly calculate times for
061             *             years <> 2000.
062             */
063            public double getUTCSunrise(AstronomicalCalendar astronomicalCalendar, double zenith, boolean adjustForElevation) {
064    //              if (astronomicalCalendar.getCalendar().get(Calendar.YEAR) == 2000) {
065    //                      throw new ZmanimException(
066    //                                      "JSuntimeCalculator can not calculate times for the year 2000. Please try a date with a different year.");
067    //              }
068    
069                    if (adjustForElevation) {
070                            zenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation()
071                                            .getElevation());
072                    } else {
073                            zenith = adjustZenith(zenith, 0);
074                    }
075                    double timeMins = morningPhenomenon(dateToJulian(astronomicalCalendar
076                                    .getCalendar()), astronomicalCalendar.getGeoLocation()
077                                    .getLatitude(), -astronomicalCalendar.getGeoLocation()
078                                    .getLongitude(), zenith);
079                    return timeMins / 60;
080            }
081    
082            /**
083             * @deprecated  This class is based on the NOAA algorithm but does not return calculations that match the NOAAA algorithm JavaScript implementation. The calculations are about 2 minutes off. This call has been replaced by the NOAACalculator class.
084             * @see net.sourceforge.zmanim.util.NOAACalculator#getUTCSunset(AstronomicalCalendar, double, boolean)
085             * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunset(AstronomicalCalendar,
086             *      double, boolean)
087             * @throws ZmanimException
088             *             if the year entered == 2000. This calculator can't properly
089             *             deal with the year 2000. It can properly calculate times for
090             *             years <> 2000.
091             */
092            public double getUTCSunset(AstronomicalCalendar astronomicalCalendar, double zenith, boolean adjustForElevation) {
093    //              if (astronomicalCalendar.getCalendar().get(Calendar.YEAR) == 2000) {
094    //                      throw new ZmanimException(
095    //                                      "JSuntimeCalculator can not calculate times for the year 2000. Please try a date with a different year.");
096    //              }
097    
098                    if (adjustForElevation) {
099                            zenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation()
100                                            .getElevation());
101                    } else {
102                            zenith = adjustZenith(zenith, 0);
103                    }
104                    double timeMins = eveningPhenomenon(dateToJulian(astronomicalCalendar
105                                    .getCalendar()), astronomicalCalendar.getGeoLocation()
106                                    .getLatitude(), -astronomicalCalendar.getGeoLocation()
107                                    .getLongitude(), zenith);
108                    return timeMins / 60;
109            }
110    
111            /**
112             * Calculate the UTC of a morning phenomenon for the given day at the given
113             * latitude and longitude on Earth
114             *
115             * @param julian
116             *            Julian day
117             * @param latitude
118             *            latitude of observer in degrees
119             * @param longitude
120             *            longitude of observer in degrees
121             * @param zenithDistance
122             *            one of Sun.SUNRISE_SUNSET_ZENITH_DISTANCE,
123             *            Sun.CIVIL_TWILIGHT_ZENITH_DISTANCE,
124             *            Sun.NAUTICAL_TWILIGHT_ZENITH_DISTANCE,
125             *            Sun.ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE.
126             * @return time in minutes from zero Z
127             */
128            private static double morningPhenomenon(double julian, double latitude,
129                            double longitude, double zenithDistance) {
130                    double t = julianDayToJulianCenturies(julian);
131                    double eqtime = equationOfTime(t);
132                    double solarDec = sunDeclination(t);
133                    double hourangle = hourAngleMorning(latitude, solarDec, zenithDistance);
134                    double delta = longitude - Math.toDegrees(hourangle);
135                    double timeDiff = 4 * delta;
136                    double timeUTC = 720 + timeDiff - eqtime;
137    
138                    // Second pass includes fractional julian day in gamma calc
139                    double newt = julianDayToJulianCenturies(julianCenturiesToJulianDay(t)
140                                    + timeUTC / 1440);
141                    eqtime = equationOfTime(newt);
142                    solarDec = sunDeclination(newt);
143                    hourangle = hourAngleMorning(latitude, solarDec, zenithDistance);
144                    delta = longitude - Math.toDegrees(hourangle);
145                    timeDiff = 4 * delta;
146    
147                    double morning = 720 + timeDiff - eqtime;
148                    return morning;
149            }
150    
151            /**
152             * Calculate the UTC of an evening phenomenon for the given day at the given
153             * latitude and longitude on Earth
154             *
155             * @param julian
156             *            Julian day
157             * @param latitude
158             *            latitude of observer in degrees
159             * @param longitude
160             *            longitude of observer in degrees
161             * @param zenithDistance
162             *            one of Sun.SUNRISE_SUNSET_ZENITH_DISTANCE,
163             *            Sun.CIVIL_TWILIGHT_ZENITH_DISTANCE,
164             *            Sun.NAUTICAL_TWILIGHT_ZENITH_DISTANCE,
165             *            Sun.ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE.
166             * @return time in minutes from zero Z
167             */
168            private static double eveningPhenomenon(double julian, double latitude,
169                            double longitude, double zenithDistance) {
170                    double t = julianDayToJulianCenturies(julian);
171    
172                    // First calculates sunrise and approx length of day
173                    double eqtime = equationOfTime(t);
174                    double solarDec = sunDeclination(t);
175                    double hourangle = hourAngleEvening(latitude, solarDec, zenithDistance);
176    
177                    double delta = longitude - Math.toDegrees(hourangle);
178                    double timeDiff = 4 * delta;
179                    double timeUTC = 720 + timeDiff - eqtime;
180    
181                    // first pass used to include fractional day in gamma calc
182                    double newt = julianDayToJulianCenturies(julianCenturiesToJulianDay(t)
183                                    + timeUTC / 1440);
184                    eqtime = equationOfTime(newt);
185                    solarDec = sunDeclination(newt);
186                    hourangle = hourAngleEvening(latitude, solarDec, zenithDistance);
187    
188                    delta = longitude - Math.toDegrees(hourangle);
189                    timeDiff = 4 * delta;
190    
191                    double evening = 720 + timeDiff - eqtime;
192                    return evening;
193            }
194    
195            private static double dateToJulian(Calendar date) {
196                    int year = date.get(Calendar.YEAR);
197                    int month = date.get(Calendar.MONTH) + 1;
198                    int day = date.get(Calendar.DAY_OF_MONTH);
199                    int hour = date.get(Calendar.HOUR_OF_DAY);
200                    int minute = date.get(Calendar.MINUTE);
201                    int second = date.get(Calendar.SECOND);
202    
203                    double extra = (100.0 * year) + month - 190002.5;
204                    double JD = (367.0 * year)
205                                    - (Math
206                                                    .floor(7.0 * (year + Math.floor((month + 9.0) / 12.0)) / 4.0))
207                                    + Math.floor((275.0 * month) / 9.0) + day
208                                    + ((hour + ((minute + (second / 60.0)) / 60.0)) / 24.0)
209                                    + 1721013.5 - ((0.5 * extra) / Math.abs(extra)) + 0.5;
210                    return JD;
211            }
212    
213            /**
214             * Convert Julian Day to centuries since J2000.0
215             *
216             * @param julian
217             *            The Julian Day to convert
218             * @return the value corresponding to the Julian Day
219             */
220            private static double julianDayToJulianCenturies(double julian) {
221                    return (julian - 2451545) / 36525;
222            }
223    
224            /**
225             * Convert centuries since J2000.0 to Julian Day
226             *
227             * @param t
228             *            Number of Julian centuries since J2000.0
229             * @return The Julian Day corresponding to the value of t
230             */
231            private static double julianCenturiesToJulianDay(double t) {
232                    return (t * 36525) + 2451545;
233            }
234    
235            /**
236             * Calculate the difference between true solar time and mean solar time
237             *
238             * @param t
239             *            Number of Julian centuries since J2000.0
240             * @return
241             */
242            private static double equationOfTime(double t) {
243                    double epsilon = obliquityCorrection(t);
244                    double l0 = geomMeanLongSun(t);
245                    double e = eccentricityOfEarthsOrbit(t);
246                    double m = geometricMeanAnomalyOfSun(t);
247                    double y = Math.pow((Math.tan(Math.toRadians(epsilon) / 2)), 2);
248    
249                    double eTime = y * Math.sin(2 * Math.toRadians(l0)) - 2 * e
250                                    * Math.sin(Math.toRadians(m)) + 4 * e * y
251                                    * Math.sin(Math.toRadians(m))
252                                    * Math.cos(2 * Math.toRadians(l0)) - 0.5 * y * y
253                                    * Math.sin(4 * Math.toRadians(l0)) - 1.25 * e * e
254                                    * Math.sin(2 * Math.toRadians(m));
255                    return Math.toDegrees(eTime) * 4;
256            }
257    
258            /**
259             * Calculate the declination of the sun
260             *
261             * @param t
262             *            Number of Julian centuries since J2000.0
263             * @return The Sun's declination in degrees
264             */
265            private static double sunDeclination(double t) {
266                    double e = obliquityCorrection(t);
267                    double lambda = sunsApparentLongitude(t);
268    
269                    double sint = Math.sin(Math.toRadians(e))
270                                    * Math.sin(Math.toRadians(lambda));
271                    return Math.toDegrees(Math.asin(sint));
272            }
273    
274            /**
275             * calculate the hour angle of the sun for a morning phenomenon for the
276             * given latitude
277             *
278             * @param lat
279             *            Latitude of the observer in degrees
280             * @param solarDec
281             *            declination of the sun in degrees
282             * @param zenithDistance
283             *            zenith distance of the sun in degrees
284             * @return hour angle of sunrise in radians
285             */
286            private static double hourAngleMorning(double lat, double solarDec,
287                            double zenithDistance) {
288                    return (Math.acos(Math.cos(Math.toRadians(zenithDistance))
289                                    / (Math.cos(Math.toRadians(lat)) * Math.cos(Math
290                                                    .toRadians(solarDec))) - Math.tan(Math.toRadians(lat))
291                                    * Math.tan(Math.toRadians(solarDec))));
292            }
293    
294            /**
295             * Calculate the hour angle of the sun for an evening phenomenon for the
296             * given latitude
297             *
298             * @param lat
299             *            Latitude of the observer in degrees
300             * @param solarDec
301             *            declination of the Sun in degrees
302             * @param zenithDistance
303             *            zenith distance of the sun in degrees
304             * @return hour angle of sunset in radians
305             */
306            private static double hourAngleEvening(double lat, double solarDec,
307                            double zenithDistance) {
308                    return -hourAngleMorning(lat, solarDec, zenithDistance);
309            }
310    
311            /**
312             * Calculate the corrected obliquity of the ecliptic
313             *
314             * @param t
315             *            Number of Julian centuries since J2000.0
316             * @return Corrected obliquity in degrees
317             */
318            private static double obliquityCorrection(double t) {
319                    return meanObliquityOfEcliptic(t) + 0.00256
320                                    * Math.cos(Math.toRadians(125.04 - 1934.136 * t));
321            }
322    
323            /**
324             * Calculate the mean obliquity of the ecliptic
325             *
326             * @param t
327             *            Number of Julian centuries since J2000.0
328             * @return Mean obliquity in degrees
329             */
330            private static double meanObliquityOfEcliptic(double t) {
331                    return 23 + (26 + (21.448 - t
332                                    * (46.815 + t * (0.00059 - t * (0.001813))) / 60)) / 60;
333            }
334    
335            /**
336             * Calculate the geometric mean longitude of the sun
337             *
338             * @param t
339             *            number of Julian centuries since J2000.0
340             * @return the geometric mean longitude of the sun in degrees
341             */
342            private static double geomMeanLongSun(double t) {
343                    double l0 = 280.46646 + t * (36000.76983 + 0.0003032 * t);
344    
345                    while ((l0 >= 0) && (l0 <= 360)) {
346                            if (l0 > 360) {
347                                    l0 = l0 - 360;
348                            }
349    
350                            if (l0 < 0) {
351                                    l0 = l0 + 360;
352                            }
353                    }
354                    return l0;
355            }
356    
357            /**
358             * Calculate the eccentricity of Earth's orbit
359             *
360             * @param t
361             *            Number of Julian centuries since J2000.0
362             * @return the eccentricity
363             */
364            private static double eccentricityOfEarthsOrbit(double t) {
365                    return 0.016708634 - t * (0.000042037 + 0.0000001267 * t);
366            }
367    
368            /**
369             * Calculate the geometric mean anomaly of the Sun
370             *
371             * @param t
372             *            Number of Julian centuries since J2000.0
373             * @return the geometric mean anomaly of the Sun in degrees
374             */
375            private static double geometricMeanAnomalyOfSun(double t) {
376                    return 357.52911 + t * (35999.05029 - 0.0001537 * t);
377            }
378    
379            /**
380             * Calculate the apparent longitude of the sun
381             *
382             * @param t
383             *            Number of Julian centuries since J2000.0
384             * @return The apparent longitude of the Sun in degrees
385             */
386            private static double sunsApparentLongitude(double t) {
387                    return sunsTrueLongitude(t) - 0.00569 - 0.00478
388                                    * Math.sin(Math.toRadians(125.04 - 1934.136 * t));
389            }
390    
391            /**
392             * Calculate the true longitude of the sun
393             *
394             * @param t
395             *            Number of Julian centuries since J2000.0
396             * @return The Sun's true longitude in degrees
397             */
398            private static double sunsTrueLongitude(double t) {
399                    return geomMeanLongSun(t) + equationOfCentreForSun(t);
400            }
401    
402            /**
403             * Calculate the equation of centre for the Sun
404             *
405             * @param centuries
406             *            Number of Julian centuries since J2000.0
407             * @return The equation of centre for the Sun in degrees
408             */
409            private static double equationOfCentreForSun(double t) {
410                    double m = geometricMeanAnomalyOfSun(t);
411    
412                    return Math.sin(Math.toRadians(m))
413                                    * (1.914602 - t * (0.004817 + 0.000014 * t))
414                                    + Math.sin(2 * Math.toRadians(m)) * (0.019993 - 0.000101 * t)
415                                    + Math.sin(3 * Math.toRadians(m)) * 0.000289;
416            }
417    }