001    /*
002     * Zmanim Java API
003     * Copyright (C) 2004-2007 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.
025     * This calculator uses the Java algorithm written by <a
026     * href="http://www.jstot.me.uk/jsuntimes/">Jonathan Stott</a> that is based on
027     * the implementation by <a href=""http://noaa.gov">NOAA - National Oceanic and
028     * Atmospheric Administration</a>'s <a href =
029     * "http://www.srrb.noaa.gov/highlights/sunrise/sunrisehtml">Surface Radiation
030     * Research Branch</a>. NOAA's <a
031     * href="http://www.srrb.noaa.gov/highlights/sunrise/solareqns.PDF">implementation</a>
032     * is based on equations from <a
033     * href="http://www.willbell.com/math/mc1.htm">Astronomical Algorithms</a> by
034     * <a href="http://en.wikipedia.org/wiki/Jean_Meeus">Jean Meeus</a>. Jonathan's
035     * implementation was released under the GPL. Added to the algorithm is an
036     * adjustment of the zenith to account for elevation.
037     *
038     * @author Jonathan Stott 2000 - 2004
039     * @author © Eliyahu Hershfeld 2004 - 2007
040     * @version 1.1
041     * @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.
042     * @see net.sourceforge.zmanim.util.NOAACalculator
043     */
044    public class JSuntimeCalculator extends AstronomicalCalculator {
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 "US National Oceanic and Atmospheric Administration Algorithm";
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                    //if(timeMins/60 < 4) {System.out.println("zenith="+zenith);}
080                    return timeMins / 60;
081            }
082    
083            /**
084             * @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.
085             * @see net.sourceforge.zmanim.util.NOAACalculator#getUTCSunset(AstronomicalCalendar, double, boolean)
086             * @see net.sourceforge.zmanim.util.AstronomicalCalculator#getUTCSunset(AstronomicalCalendar,
087             *      double, boolean)
088             * @throws ZmanimException
089             *             if the year entered == 2000. This calculator can't properly
090             *             deal with the year 2000. It can properly calculate times for
091             *             years <> 2000.
092             */
093            public double getUTCSunset(AstronomicalCalendar astronomicalCalendar, double zenith, boolean adjustForElevation) {
094                    if (astronomicalCalendar.getCalendar().get(Calendar.YEAR) == 2000) {
095                            throw new ZmanimException(
096                                            "JSuntimeCalculator can not calculate times for the year 2000. Please try a date with a different year.");
097                    }
098    
099                    if (adjustForElevation) {
100                            zenith = adjustZenith(zenith, astronomicalCalendar.getGeoLocation()
101                                            .getElevation());
102                    } else {
103                            zenith = adjustZenith(zenith, 0);
104                    }
105                    double timeMins = eveningPhenomenon(dateToJulian(astronomicalCalendar
106                                    .getCalendar()), astronomicalCalendar.getGeoLocation()
107                                    .getLatitude(), -astronomicalCalendar.getGeoLocation()
108                                    .getLongitude(), zenith);
109                    //if(timeMins/60 > 28) {System.out.println("zenith="+zenith);}
110                    return timeMins / 60;
111            }
112    
113            /**
114             * Calculate the UTC of a morning phenomenon for the given day at the given
115             * latitude and longitude on Earth
116             *
117             * @param julian
118             *            Julian day
119             * @param latitude
120             *            latitude of observer in degrees
121             * @param longitude
122             *            longitude of observer in degrees
123             * @param zenithDistance
124             *            one of Sun.SUNRISE_SUNSET_ZENITH_DISTANCE,
125             *            Sun.CIVIL_TWILIGHT_ZENITH_DISTANCE,
126             *            Sun.NAUTICAL_TWILIGHT_ZENITH_DISTANCE,
127             *            Sun.ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE.
128             * @return time in minutes from zero Z
129             */
130            private static double morningPhenomenon(double julian, double latitude,
131                            double longitude, double zenithDistance) {
132                    double t = julianDayToJulianCenturies(julian);
133                    double eqtime = equationOfTime(t);
134                    double solarDec = sunDeclination(t);
135                    double hourangle = hourAngleMorning(latitude, solarDec, zenithDistance);
136                    double delta = longitude - Math.toDegrees(hourangle);
137                    double timeDiff = 4 * delta;
138                    double timeUTC = 720 + timeDiff - eqtime;
139    
140                    // Second pass includes fractional julian day in gamma calc
141                    double newt = julianDayToJulianCenturies(julianCenturiesToJulianDay(t)
142                                    + timeUTC / 1440);
143                    eqtime = equationOfTime(newt);
144                    solarDec = sunDeclination(newt);
145                    hourangle = hourAngleMorning(latitude, solarDec, zenithDistance);
146                    delta = longitude - Math.toDegrees(hourangle);
147                    timeDiff = 4 * delta;
148    
149                    return 720 + timeDiff - eqtime;
150            }
151    
152            /**
153             * Calculate the UTC of an evening phenomenon for the given day at the given
154             * latitude and longitude on Earth
155             *
156             * @param julian
157             *            Julian day
158             * @param latitude
159             *            latitude of observer in degrees
160             * @param longitude
161             *            longitude of observer in degrees
162             * @param zenithDistance
163             *            one of Sun.SUNRISE_SUNSET_ZENITH_DISTANCE,
164             *            Sun.CIVIL_TWILIGHT_ZENITH_DISTANCE,
165             *            Sun.NAUTICAL_TWILIGHT_ZENITH_DISTANCE,
166             *            Sun.ASTRONOMICAL_TWILIGHT_ZENITH_DISTANCE.
167             * @return time in minutes from zero Z
168             */
169            private static double eveningPhenomenon(double julian, double latitude,
170                            double longitude, double zenithDistance) {
171                    double t = julianDayToJulianCenturies(julian);
172    
173                    // First calculates sunrise and approx length of day
174                    double eqtime = equationOfTime(t);
175                    double solarDec = sunDeclination(t);
176                    double hourangle = hourAngleEvening(latitude, solarDec, zenithDistance);
177    
178                    double delta = longitude - Math.toDegrees(hourangle);
179                    double timeDiff = 4 * delta;
180                    double timeUTC = 720 + timeDiff - eqtime;
181    
182                    // first pass used to include fractional day in gamma calc
183                    double newt = julianDayToJulianCenturies(julianCenturiesToJulianDay(t)
184                                    + timeUTC / 1440);
185                    eqtime = equationOfTime(newt);
186                    solarDec = sunDeclination(newt);
187                    hourangle = hourAngleEvening(latitude, solarDec, zenithDistance);
188    
189                    delta = longitude - Math.toDegrees(hourangle);
190                    timeDiff = 4 * delta;
191    
192                    return 720 + timeDiff - eqtime;
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                    return (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            }
211    
212            /**
213             * Convert Julian Day to centuries since J2000.0
214             *
215             * @param julian
216             *            The Julian Day to convert
217             * @return the value corresponding to the Julian Day
218             */
219            private static double julianDayToJulianCenturies(double julian) {
220                    return (julian - 2451545) / 36525;
221            }
222    
223            /**
224             * Convert centuries since J2000.0 to Julian Day
225             *
226             * @param t
227             *            Number of Julian centuries since J2000.0
228             * @return The Julian Day corresponding to the value of t
229             */
230            private static double julianCenturiesToJulianDay(double t) {
231                    return (t * 36525) + 2451545;
232            }
233    
234            /**
235             * Calculate the difference between true solar time and mean solar time
236             *
237             * @param t
238             *            Number of Julian centuries since J2000.0
239             * @return
240             */
241            private static double equationOfTime(double t) {
242                    double epsilon = obliquityCorrection(t);
243                    double l0 = geomMeanLongSun(t);
244                    double e = eccentricityOfEarthsOrbit(t);
245                    double m = geometricMeanAnomalyOfSun(t);
246                    double y = Math.pow((Math.tan(Math.toRadians(epsilon) / 2)), 2);
247    
248                    double eTime = y * Math.sin(2 * Math.toRadians(l0)) - 2 * e
249                                    * Math.sin(Math.toRadians(m)) + 4 * e * y
250                                    * Math.sin(Math.toRadians(m))
251                                    * Math.cos(2 * Math.toRadians(l0)) - 0.5 * y * y
252                                    * Math.sin(4 * Math.toRadians(l0)) - 1.25 * e * e
253                                    * Math.sin(2 * Math.toRadians(m));
254    
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    
355                    return l0;
356            }
357    
358            /**
359             * Calculate the eccentricity of Earth's orbit
360             *
361             * @param t
362             *            Number of Julian centuries since J2000.0
363             * @return the eccentricity
364             */
365            private static double eccentricityOfEarthsOrbit(double t) {
366                    return 0.016708634 - t * (0.000042037 + 0.0000001267 * t);
367            }
368    
369            /**
370             * Calculate the geometric mean anomaly of the Sun
371             *
372             * @param t
373             *            Number of Julian centuries since J2000.0
374             * @return the geometric mean anomaly of the Sun in degrees
375             */
376            private static double geometricMeanAnomalyOfSun(double t) {
377                    return 357.52911 + t * (35999.05029 - 0.0001537 * t);
378            }
379    
380            /**
381             * Calculate the apparent longitude of the sun
382             *
383             * @param t
384             *            Number of Julian centuries since J2000.0
385             * @return The apparent longitude of the Sun in degrees
386             */
387            private static double sunsApparentLongitude(double t) {
388                    return sunsTrueLongitude(t) - 0.00569 - 0.00478
389                                    * Math.sin(Math.toRadians(125.04 - 1934.136 * t));
390            }
391    
392            /**
393             * Calculate the true longitude of the sun
394             *
395             * @param t
396             *            Number of Julian centuries since J2000.0
397             * @return The Sun's true longitude in degrees
398             */
399            private static double sunsTrueLongitude(double t) {
400                    return geomMeanLongSun(t) + equationOfCentreForSun(t);
401            }
402    
403            /**
404             * Calculate the equation of centre for the Sun
405             *
406             * @param centuries
407             *            Number of Julian centuries since J2000.0
408             * @return The equation of centre for the Sun in degrees
409             */
410            private static double equationOfCentreForSun(double t) {
411                    double m = geometricMeanAnomalyOfSun(t);
412    
413                    return Math.sin(Math.toRadians(m))
414                                    * (1.914602 - t * (0.004817 + 0.000014 * t))
415                                    + Math.sin(2 * Math.toRadians(m)) * (0.019993 - 0.000101 * t)
416                                    + Math.sin(3 * Math.toRadians(m)) * 0.000289;
417            }
418    }