001/*
002 * Zmanim Java API
003 * Copyright (C) 2004-2023 Eliyahu Hershfeld
004 *
005 * This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General
006 * Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option)
007 * any later version.
008 *
009 * This library is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY; without even the implied
010 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
011 * details.
012 * You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to
013 * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA,
014 * or connect to: https://www.gnu.org/licenses/old-licenses/lgpl-2.1.html
015 */
016package com.kosherjava.zmanim.util;
017
018import java.util.Calendar;
019
020/**
021 * Implementation of sunrise and sunset methods to calculate astronomical times based on the <a
022 * href="https://noaa.gov">NOAA</a> algorithm. This calculator uses the Java algorithm based on the implementation by <a
023 * href="https://noaa.gov">NOAA - National Oceanic and Atmospheric Administration</a>'s <a href =
024 * "https://www.srrb.noaa.gov/highlights/sunrise/sunrise.html">Surface Radiation Research Branch</a>. NOAA's <a
025 * href="https://www.srrb.noaa.gov/highlights/sunrise/solareqns.PDF">implementation</a> is based on equations from <a
026 * href="https://www.amazon.com/Astronomical-Table-Sun-Moon-Planets/dp/1942675038/">Astronomical Algorithms</a> by <a
027 * href="https://en.wikipedia.org/wiki/Jean_Meeus">Jean Meeus</a>. Added to the algorithm is an adjustment of the zenith
028 * to account for elevation. The algorithm can be found in the <a
029 * href="https://en.wikipedia.org/wiki/Sunrise_equation">Wikipedia Sunrise Equation</a> article.
030 * 
031 * @author &copy; Eliyahu Hershfeld 2011 - 2023
032 */
033public class NOAACalculator extends AstronomicalCalculator {
034        /**
035         * The <a href="https://en.wikipedia.org/wiki/Julian_day">Julian day</a> of January 1, 2000, known as
036         * <a href="https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
037         */
038        private static final double JULIAN_DAY_JAN_1_2000 = 2451545.0;
039
040        /**
041         * Julian days per century.
042         */
043        private static final double JULIAN_DAYS_PER_CENTURY = 36525.0;
044
045        /**
046         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getCalculatorName()
047         */
048        public String getCalculatorName() {
049                return "US National Oceanic and Atmospheric Administration Algorithm";
050        }
051
052        /**
053         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCSunrise(Calendar, GeoLocation, double, boolean)
054         */
055        public double getUTCSunrise(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) {
056                double elevation = adjustForElevation ? geoLocation.getElevation() : 0;
057                double adjustedZenith = adjustZenith(zenith, elevation);
058
059                double sunrise = getSunriseUTC(getJulianDay(calendar), geoLocation.getLatitude(), -geoLocation.getLongitude(),
060                                adjustedZenith);
061                sunrise = sunrise / 60;
062
063                // ensure that the time is >= 0 and < 24
064                while (sunrise < 0.0) {
065                        sunrise += 24.0;
066                }
067                while (sunrise >= 24.0) {
068                        sunrise -= 24.0;
069                }
070                return sunrise;
071        }
072
073        /**
074         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCSunset(Calendar, GeoLocation, double, boolean)
075         */
076        public double getUTCSunset(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) {
077                double elevation = adjustForElevation ? geoLocation.getElevation() : 0;
078                double adjustedZenith = adjustZenith(zenith, elevation);
079
080                double sunset = getSunsetUTC(getJulianDay(calendar), geoLocation.getLatitude(), -geoLocation.getLongitude(),
081                                adjustedZenith);
082                sunset = sunset / 60;
083
084                // ensure that the time is >= 0 and < 24
085                while (sunset < 0.0) {
086                        sunset += 24.0;
087                }
088                while (sunset >= 24.0) {
089                        sunset -= 24.0;
090                }
091                return sunset;
092        }
093
094        /**
095         * Return the <a href="https://en.wikipedia.org/wiki/Julian_day">Julian day</a> from a Java Calendar
096         * 
097         * @param calendar
098         *            The Java Calendar
099         * @return the Julian day corresponding to the date Note: Number is returned for start of day. Fractional days
100         *         should be added later.
101         */
102        private static double getJulianDay(Calendar calendar) {
103                int year = calendar.get(Calendar.YEAR);
104                int month = calendar.get(Calendar.MONTH) + 1;
105                int day = calendar.get(Calendar.DAY_OF_MONTH);
106                if (month <= 2) {
107                        year -= 1;
108                        month += 12;
109                }
110                int a = year / 100;
111                int b = 2 - a + a / 4;
112
113                return Math.floor(365.25 * (year + 4716)) + Math.floor(30.6001 * (month + 1)) + day + b - 1524.5;
114        }
115
116        /**
117         * Convert <a href="https://en.wikipedia.org/wiki/Julian_day">Julian day</a> to centuries since <a href=
118         * "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
119         * 
120         * @param julianDay
121         *            the Julian Day to convert
122         * @return the centuries since 2000 Julian corresponding to the Julian Day
123         */
124        private static double getJulianCenturiesFromJulianDay(double julianDay) {
125                return (julianDay - JULIAN_DAY_JAN_1_2000) / JULIAN_DAYS_PER_CENTURY;
126        }
127
128        /**
129         * Convert centuries since <a href="https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a> to
130         * <a href="https://en.wikipedia.org/wiki/Julian_day">Julian day</a>.
131         * 
132         * @param julianCenturies
133         *            the number of Julian centuries since <a href=
134         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
135         * @return the Julian Day corresponding to the Julian centuries passed in
136         */
137        private static double getJulianDayFromJulianCenturies(double julianCenturies) {
138                return julianCenturies * JULIAN_DAYS_PER_CENTURY + JULIAN_DAY_JAN_1_2000;
139        }
140
141        /**
142         * Returns the Geometric <a href="https://en.wikipedia.org/wiki/Mean_longitude">Mean Longitude</a> of the Sun.
143         * 
144         * @param julianCenturies
145         *            the number of Julian centuries since <a href=
146         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
147         * @return the Geometric Mean Longitude of the Sun in degrees
148         */
149        private static double getSunGeometricMeanLongitude(double julianCenturies) {
150                double longitude = 280.46646 + julianCenturies * (36000.76983 + 0.0003032 * julianCenturies);
151                while (longitude > 360.0) {
152                        longitude -= 360.0;
153                }
154                while (longitude < 0.0) {
155                        longitude += 360.0;
156                }
157
158                return longitude; // in degrees
159        }
160
161        /**
162         * Returns the Geometric <a href="https://en.wikipedia.org/wiki/Mean_anomaly">Mean Anomaly</a> of the Sun.
163         * 
164         * @param julianCenturies
165         *            the number of Julian centuries since <a href=
166         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
167         * @return the Geometric Mean Anomaly of the Sun in degrees
168         */
169        private static double getSunGeometricMeanAnomaly(double julianCenturies) {
170                return 357.52911 + julianCenturies * (35999.05029 - 0.0001537 * julianCenturies); // in degrees
171        }
172
173        /**
174         * Return the <a href="https://en.wikipedia.org/wiki/Eccentricity_%28orbit%29">eccentricity of earth's orbit</a>.
175         * 
176         * @param julianCenturies
177         *            the number of Julian centuries since <a href=
178         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
179         * @return the unitless eccentricity
180         */
181        private static double getEarthOrbitEccentricity(double julianCenturies) {
182                return 0.016708634 - julianCenturies * (0.000042037 + 0.0000001267 * julianCenturies); // unitless
183        }
184
185        /**
186         * Returns the <a href="https://en.wikipedia.org/wiki/Equation_of_the_center">equation of center</a> for the sun.
187         * 
188         * @param julianCenturies
189         *            the number of Julian centuries since <a href=
190         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
191         * @return the equation of center for the sun in degrees
192         */
193        private static double getSunEquationOfCenter(double julianCenturies) {
194                double m = getSunGeometricMeanAnomaly(julianCenturies);
195
196                double mrad = Math.toRadians(m);
197                double sinm = Math.sin(mrad);
198                double sin2m = Math.sin(mrad + mrad);
199                double sin3m = Math.sin(mrad + mrad + mrad);
200
201                return sinm * (1.914602 - julianCenturies * (0.004817 + 0.000014 * julianCenturies)) + sin2m
202                                * (0.019993 - 0.000101 * julianCenturies) + sin3m * 0.000289;// in degrees
203        }
204
205        /**
206         * Return the <a href="https://en.wikipedia.org/wiki/True_longitude">true longitude</a> of the sun.
207         * 
208         * @param julianCenturies
209         *            the number of Julian centuries since <a href=
210         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
211         * @return the sun's true longitude in degrees
212         */
213        private static double getSunTrueLongitude(double julianCenturies) {
214                double sunLongitude = getSunGeometricMeanLongitude(julianCenturies);
215                double center = getSunEquationOfCenter(julianCenturies);
216
217                return sunLongitude + center; // in degrees
218        }
219
220        // /**
221        // * Returns the <a href="https://en.wikipedia.org/wiki/True_anomaly">true anamoly</a> of the sun.
222        // *
223        // * @param julianCenturies
224        // * the number of Julian centuries since J2000.0
225        // * @return the sun's true anamoly in degrees
226        // */
227        // private static double getSunTrueAnomaly(double julianCenturies) {
228        // double meanAnomaly = getSunGeometricMeanAnomaly(julianCenturies);
229        // double equationOfCenter = getSunEquationOfCenter(julianCenturies);
230        //
231        // return meanAnomaly + equationOfCenter; // in degrees
232        // }
233
234        /**
235         * Return the <a href="https://en.wikipedia.org/wiki/Apparent_longitude">apparent longitude</a> of the sun.
236         * 
237         * @param julianCenturies
238         *            the number of Julian centuries since <a href=
239         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
240         * @return sun's apparent longitude in degrees
241         */
242        private static double getSunApparentLongitude(double julianCenturies) {
243                double sunTrueLongitude = getSunTrueLongitude(julianCenturies);
244
245                double omega = 125.04 - 1934.136 * julianCenturies;
246                double lambda = sunTrueLongitude - 0.00569 - 0.00478 * Math.sin(Math.toRadians(omega));
247                return lambda; // in degrees
248        }
249
250        /**
251         * Returns the mean <a href="https://en.wikipedia.org/wiki/Axial_tilt">obliquity of the ecliptic</a> (Axial tilt).
252         * 
253         * @param julianCenturies
254         *            the number of Julian centuries since <a href=
255         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
256         * @return the mean obliquity in degrees
257         */
258        private static double getMeanObliquityOfEcliptic(double julianCenturies) {
259                double seconds = 21.448 - julianCenturies
260                                * (46.8150 + julianCenturies * (0.00059 - julianCenturies * (0.001813)));
261                return 23.0 + (26.0 + (seconds / 60.0)) / 60.0; // in degrees
262        }
263
264        /**
265         * Returns the corrected <a href="https://en.wikipedia.org/wiki/Axial_tilt">obliquity of the ecliptic</a> (Axial
266         * tilt).
267         * 
268         * @param julianCenturies
269         *            the number of Julian centuries since <a href=
270         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
271         * @return the corrected obliquity in degrees
272         */
273        private static double getObliquityCorrection(double julianCenturies) {
274                double obliquityOfEcliptic = getMeanObliquityOfEcliptic(julianCenturies);
275
276                double omega = 125.04 - 1934.136 * julianCenturies;
277                return obliquityOfEcliptic + 0.00256 * Math.cos(Math.toRadians(omega)); // in degrees
278        }
279
280        /**
281         * Return the <a href="https://en.wikipedia.org/wiki/Declination">declination</a> of the sun.
282         * 
283         * @param julianCenturies
284         *            the number of Julian centuries since <a href=
285         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
286         * @return
287         *            the sun's declination in degrees
288         */
289        private static double getSunDeclination(double julianCenturies) {
290                double obliquityCorrection = getObliquityCorrection(julianCenturies);
291                double lambda = getSunApparentLongitude(julianCenturies);
292
293                double sint = Math.sin(Math.toRadians(obliquityCorrection)) * Math.sin(Math.toRadians(lambda));
294                double theta = Math.toDegrees(Math.asin(sint));
295                return theta; // in degrees
296        }
297
298        /**
299         * Return the <a href="https://en.wikipedia.org/wiki/Equation_of_time">Equation of Time</a> - the difference between
300         * true solar time and mean solar time
301         * 
302         * @param julianCenturies
303         *            the number of Julian centuries since <a href=
304         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
305         * @return equation of time in minutes of time
306         */
307        private static double getEquationOfTime(double julianCenturies) {
308                double epsilon = getObliquityCorrection(julianCenturies);
309                double geomMeanLongSun = getSunGeometricMeanLongitude(julianCenturies);
310                double eccentricityEarthOrbit = getEarthOrbitEccentricity(julianCenturies);
311                double geomMeanAnomalySun = getSunGeometricMeanAnomaly(julianCenturies);
312
313                double y = Math.tan(Math.toRadians(epsilon) / 2.0);
314                y *= y;
315
316                double sin2l0 = Math.sin(2.0 * Math.toRadians(geomMeanLongSun));
317                double sinm = Math.sin(Math.toRadians(geomMeanAnomalySun));
318                double cos2l0 = Math.cos(2.0 * Math.toRadians(geomMeanLongSun));
319                double sin4l0 = Math.sin(4.0 * Math.toRadians(geomMeanLongSun));
320                double sin2m = Math.sin(2.0 * Math.toRadians(geomMeanAnomalySun));
321
322                double equationOfTime = y * sin2l0 - 2.0 * eccentricityEarthOrbit * sinm + 4.0 * eccentricityEarthOrbit * y
323                                * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * eccentricityEarthOrbit * eccentricityEarthOrbit * sin2m;
324                return Math.toDegrees(equationOfTime) * 4.0; // in minutes of time
325        }
326
327        /**
328         * Return the <a href="https://en.wikipedia.org/wiki/Hour_angle">hour angle</a> of the sun in
329         * <a href="https://en.wikipedia.org/wiki/Radian">radians</a> 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         * @param zenith
336         *            the zenith
337         * @return hour angle of sunrise in <a href="https://en.wikipedia.org/wiki/Radian">radians</a>
338         */
339        private static double getSunHourAngleAtSunrise(double lat, double solarDec, double zenith) {
340                double latRad = Math.toRadians(lat);
341                double sdRad = Math.toRadians(solarDec);
342
343                return (Math.acos(Math.cos(Math.toRadians(zenith)) / (Math.cos(latRad) * Math.cos(sdRad)) - Math.tan(latRad)
344                                * Math.tan(sdRad))); // in radians
345        }
346
347        /**
348         * Returns the <a href="https://en.wikipedia.org/wiki/Hour_angle">hour angle</a> of the sun in <a href=
349         * "https://en.wikipedia.org/wiki/Radian">radians</a>at sunset for the latitude.
350         * @todo use - {@link #getSunHourAngleAtSunrise(double, double, double)} implementation to avoid
351         * duplication of code.
352         * 
353         * @param lat
354         *            the latitude of observer in degrees
355         * @param solarDec
356         *            the declination angle of sun in degrees
357         * @param zenith
358         *            the zenith
359         * @return the hour angle of sunset in <a href="https://en.wikipedia.org/wiki/Radian">radians</a>
360         */
361        private static double getSunHourAngleAtSunset(double lat, double solarDec, double zenith) {
362                double latRad = Math.toRadians(lat);
363                double sdRad = Math.toRadians(solarDec);
364
365                double hourAngle = (Math.acos(Math.cos(Math.toRadians(zenith)) / (Math.cos(latRad) * Math.cos(sdRad))
366                                - Math.tan(latRad) * Math.tan(sdRad)));
367                return -hourAngle; // in radians
368        }
369
370        /**
371         * Return the <a href="https://en.wikipedia.org/wiki/Celestial_coordinate_system">Solar Elevation</a> for the
372         * horizontal coordinate system at the given location at the given time. Can be negative if the sun is below the
373         * horizon. Not corrected for altitude.
374         * 
375         * @param cal
376         *            time of calculation
377         * @param lat
378         *            latitude of location for calculation
379         * @param lon
380         *            longitude of location for calculation
381         * @return solar elevation in degrees - horizon is 0 degrees, civil twilight is -6 degrees
382         */
383
384        public static double getSolarElevation(Calendar cal, double lat, double lon) {
385                double julianDay = getJulianDay(cal);
386                double julianCenturies = getJulianCenturiesFromJulianDay(julianDay);
387
388                Double eot = getEquationOfTime(julianCenturies);
389
390                double longitude = (cal.get(Calendar.HOUR_OF_DAY) + 12.0)
391                                + (cal.get(Calendar.MINUTE) + eot + cal.get(Calendar.SECOND) / 60.0) / 60.0;
392
393                longitude = -(longitude * 360.0 / 24.0) % 360.0;
394                double hourAngle_rad = Math.toRadians(lon - longitude);
395                double declination = getSunDeclination(julianCenturies);
396                double dec_rad = Math.toRadians(declination);
397                double lat_rad = Math.toRadians(lat);
398                return Math.toDegrees(Math.asin((Math.sin(lat_rad) * Math.sin(dec_rad))
399                                + (Math.cos(lat_rad) * Math.cos(dec_rad) * Math.cos(hourAngle_rad))));
400
401        }
402
403        /**
404         * Return the <a href="https://en.wikipedia.org/wiki/Celestial_coordinate_system">Solar Azimuth</a> for the
405         * horizontal coordinate system at the given location at the given time. Not corrected for altitude. True south is 0
406         * degrees.
407         * 
408         * @param cal
409         *            time of calculation
410         * @param lat
411         *            latitude of location for calculation
412         * @param lon
413         *            longitude of location for calculation
414         * @return the solar azimuth
415         */
416
417        public static double getSolarAzimuth(Calendar cal, double lat, double lon) {
418                double julianDay = getJulianDay(cal);
419                double julianCenturies = getJulianCenturiesFromJulianDay(julianDay);
420
421                Double eot = getEquationOfTime(julianCenturies);
422
423                double longitude = (cal.get(Calendar.HOUR_OF_DAY) + 12.0)
424                                + (cal.get(Calendar.MINUTE) + eot + cal.get(Calendar.SECOND) / 60.0) / 60.0;
425
426                longitude = -(longitude * 360.0 / 24.0) % 360.0;
427                double hourAngle_rad = Math.toRadians(lon - longitude);
428                double declination = getSunDeclination(julianCenturies);
429                double dec_rad = Math.toRadians(declination);
430                double lat_rad = Math.toRadians(lat);
431
432                return Math.toDegrees(Math.atan(Math.sin(hourAngle_rad)
433                                / ((Math.cos(hourAngle_rad) * Math.sin(lat_rad)) - (Math.tan(dec_rad) * Math.cos(lat_rad)))))+180;
434
435        }
436
437        /**
438         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC)
439         * of sunrise for the given day at the given location on earth.
440         * 
441         * @param julianDay
442         *            the Julian day
443         * @param latitude
444         *            the latitude of observer in degrees
445         * @param longitude
446         *            the longitude of observer in degrees
447         * @param zenith
448         *            the zenith
449         * @return the time in minutes from zero UTC
450         */
451        private static double getSunriseUTC(double julianDay, double latitude, double longitude, double zenith) {
452                double julianCenturies = getJulianCenturiesFromJulianDay(julianDay);
453
454                // Find the time of solar noon at the location, and use that declination. This is better than start of the
455                // Julian day
456
457                double noonmin = getSolarNoonUTC(julianCenturies, longitude);
458                double tnoon = getJulianCenturiesFromJulianDay(julianDay + noonmin / 1440.0);
459
460                // First pass to approximate sunrise (using solar noon)
461
462                double eqTime = getEquationOfTime(tnoon);
463                double solarDec = getSunDeclination(tnoon);
464                double hourAngle = getSunHourAngleAtSunrise(latitude, solarDec, zenith);
465
466                double delta = longitude - Math.toDegrees(hourAngle);
467                double timeDiff = 4 * delta; // in minutes of time
468                double timeUTC = 720 + timeDiff - eqTime; // in minutes
469
470                // Second pass includes fractional Julian Day in gamma calc
471
472                double newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + timeUTC
473                                / 1440.0);
474                eqTime = getEquationOfTime(newt);
475                solarDec = getSunDeclination(newt);
476                hourAngle = getSunHourAngleAtSunrise(latitude, solarDec, zenith);
477                delta = longitude - Math.toDegrees(hourAngle);
478                timeDiff = 4 * delta;
479                timeUTC = 720 + timeDiff - eqTime; // in minutes
480                return timeUTC;
481        }
482        
483        /**
484         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC)
485         * of <a href="https://en.wikipedia.org/wiki/Noon#Solar_noon">solar noon</a> for the given day at the given location
486         * on earth. This implementation returns true solar noon as opposed to the time halfway between sunrise and sunset.
487         * Other calculators may return a more simplified calculation of halfway between sunrise and sunset. See <a href=
488         * "https://kosherjava.com/2020/07/02/definition-of-chatzos/">The Definition of <em>Chatzos</em></a> for details on
489         * solar noon calculations.
490         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCNoon(Calendar, GeoLocation)
491         * @see #getSolarNoonUTC(double, double)
492         * 
493         * @param calendar
494         *            The Calendar representing the date to calculate solar noon for
495         * @param geoLocation
496         *            The location information used for astronomical calculating sun times. This class uses only requires
497         *            the longitude for calculating noon since it is the same time anywhere along the longitude line.
498         * @return the time in minutes from zero UTC
499         */
500        public double getUTCNoon(Calendar calendar, GeoLocation geoLocation) {
501                double julianDay = getJulianDay(calendar);
502                double julianCenturies = getJulianCenturiesFromJulianDay(julianDay);
503                
504                double noon = getSolarNoonUTC(julianCenturies, -geoLocation.getLongitude());
505                noon = noon / 60;
506
507                // ensure that the time is >= 0 and < 24
508                while (noon < 0.0) {
509                        noon += 24.0;
510                }
511                while (noon >= 24.0) {
512                        noon -= 24.0;
513                }
514                return noon;
515        }
516
517        /**
518         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC)
519         * of of <a href="http://en.wikipedia.org/wiki/Noon#Solar_noon">solar noon</a> for the given day at the given location
520         * on earth.
521         * 
522         * @param julianCenturies
523         *            the number of Julian centuries since <a href=
524         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
525         * @param longitude
526         *            the longitude of observer in degrees
527         * 
528         * @return the time in minutes from zero UTC
529         * 
530         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCNoon(Calendar, GeoLocation)
531         * @see #getUTCNoon(Calendar, GeoLocation)
532         */
533        private static double getSolarNoonUTC(double julianCenturies, double longitude) {
534                // First pass uses approximate solar noon to calculate equation of time
535                double tnoon = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + longitude
536                                / 360.0);
537                double eqTime = getEquationOfTime(tnoon);
538                double solNoonUTC = 720 + (longitude * 4) - eqTime; // min
539
540                double newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) - 0.5
541                                + solNoonUTC / 1440.0);
542
543                eqTime = getEquationOfTime(newt);
544                return 720 + (longitude * 4) - eqTime; // min
545        }
546
547        /**
548         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC)
549         * of sunset for the given day at the given location on earth.
550         * 
551         * @param julianDay
552         *            the Julian day
553         * @param latitude
554         *            the latitude of observer in degrees
555         * @param longitude
556         *            longitude of observer in degrees
557         * @param zenith
558         *            zenith
559         * @return the time in minutes from zero Universal Coordinated Time (UTC)
560         */
561        private static double getSunsetUTC(double julianDay, double latitude, double longitude, double zenith) {
562                double julianCenturies = getJulianCenturiesFromJulianDay(julianDay);
563
564                // Find the time of solar noon at the location, and use that declination. This is better than start of the
565                // Julian day
566
567                double noonmin = getSolarNoonUTC(julianCenturies, longitude);
568                double tnoon = getJulianCenturiesFromJulianDay(julianDay + noonmin / 1440.0);
569
570                // First calculates sunrise and approx length of day
571
572                double eqTime = getEquationOfTime(tnoon);
573                double solarDec = getSunDeclination(tnoon);
574                double hourAngle = getSunHourAngleAtSunset(latitude, solarDec, zenith);
575
576                double delta = longitude - Math.toDegrees(hourAngle);
577                double timeDiff = 4 * delta;
578                double timeUTC = 720 + timeDiff - eqTime;
579
580                // Second pass includes fractional Julian Day in gamma calc
581
582                double newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + timeUTC
583                                / 1440.0);
584                eqTime = getEquationOfTime(newt);
585                solarDec = getSunDeclination(newt);
586                hourAngle = getSunHourAngleAtSunset(latitude, solarDec, zenith);
587
588                delta = longitude - Math.toDegrees(hourAngle);
589                timeDiff = 4 * delta;
590                timeUTC = 720 + timeDiff - eqTime; // in minutes
591                return timeUTC;
592        }
593}