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