001/*
002 * Zmanim Java API
003 * Copyright (C) 2004-2025 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 - 2025
032 */
033public class NOAACalculator extends AstronomicalCalculator {
034        
035        /**
036         * The <a href="https://en.wikipedia.org/wiki/Julian_day">Julian day</a> of January 1, 2000, known as
037         * <a href="https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
038         */
039        private static final double JULIAN_DAY_JAN_1_2000 = 2451545.0;
040
041        /**
042         * Julian days per century.
043         */
044        private static final double JULIAN_DAYS_PER_CENTURY = 36525.0;
045        
046        /**
047         * An <code>enum</code> to indicate what type of solar event ({@link #SUNRISE SUNRISE}, {@link #SUNSET SUNSET},
048         * {@link #NOON NOON} or {@link #MIDNIGHT MIDNIGHT}) is being calculated.
049         * .
050         */
051        protected enum SolarEvent {
052                /**SUNRISE A solar event related to sunrise*/SUNRISE, /**SUNSET A solar event related to sunset*/SUNSET,
053                /**NOON A solar event related to noon*/NOON, /**MIDNIGHT A solar event related to midnight*/MIDNIGHT
054        }
055        
056        /**
057         * Default constructor of the NOAACalculator.
058         */
059        public NOAACalculator() {
060                super();
061        }
062
063        /**
064         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getCalculatorName()
065         */
066        public String getCalculatorName() {
067                return "US National Oceanic and Atmospheric Administration Algorithm"; // Implementation of the Jean Meeus algorithm
068        }
069
070        /**
071         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCSunrise(Calendar, GeoLocation, double, boolean)
072         */
073        public double getUTCSunrise(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) {
074                double elevation = adjustForElevation ? geoLocation.getElevation() : 0;
075                double adjustedZenith = adjustZenith(zenith, elevation);
076                double sunrise = getSunRiseSetUTC(calendar, geoLocation.getLatitude(), -geoLocation.getLongitude(),
077                                adjustedZenith, SolarEvent.SUNRISE);
078                sunrise = sunrise / 60;
079                return sunrise > 0  ? sunrise % 24 : sunrise % 24 + 24; // ensure that the time is >= 0 and < 24
080        }
081
082        /**
083         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCSunset(Calendar, GeoLocation, double, boolean)
084         */
085        public double getUTCSunset(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) {
086                double elevation = adjustForElevation ? geoLocation.getElevation() : 0;
087                double adjustedZenith = adjustZenith(zenith, elevation);
088                double sunset = getSunRiseSetUTC(calendar, geoLocation.getLatitude(), -geoLocation.getLongitude(),
089                                adjustedZenith, SolarEvent.SUNSET);
090                sunset = sunset / 60;
091                return sunset > 0  ? sunset % 24 : sunset % 24 + 24; // ensure that the time is >= 0 and < 24
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                return Math.floor(365.25 * (year + 4716)) + Math.floor(30.6001 * (month + 1)) + day + b - 1524.5;
113        }
114
115        /**
116         * Convert <a href="https://en.wikipedia.org/wiki/Julian_day">Julian day</a> to centuries since <a href=
117         * "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
118         * 
119         * @param julianDay
120         *            the Julian Day to convert
121         * @return the centuries since 2000 Julian corresponding to the Julian Day
122         */
123        private static double getJulianCenturiesFromJulianDay(double julianDay) {
124                return (julianDay - JULIAN_DAY_JAN_1_2000) / JULIAN_DAYS_PER_CENTURY;
125        }
126
127        /**
128         * Returns the Geometric <a href="https://en.wikipedia.org/wiki/Mean_longitude">Mean Longitude</a> of the Sun.
129         * 
130         * @param julianCenturies
131         *            the number of Julian centuries since <a href=
132         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
133         * @return the Geometric Mean Longitude of the Sun in degrees
134         */
135        private static double getSunGeometricMeanLongitude(double julianCenturies) {
136                double longitude = 280.46646 + julianCenturies * (36000.76983 + 0.0003032 * julianCenturies);
137                return longitude > 0  ? longitude % 360 : longitude % 360 + 360; // ensure that the longitude is >= 0 and < 360
138        }
139
140        /**
141         * Returns the Geometric <a href="https://en.wikipedia.org/wiki/Mean_anomaly">Mean Anomaly</a> of the Sun in degrees.
142         * 
143         * @param julianCenturies
144         *            the number of Julian centuries since <a href=
145         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
146         * @return the Geometric Mean Anomaly of the Sun in degrees
147         */
148        private static double getSunGeometricMeanAnomaly(double julianCenturies) {
149                return 357.52911 + julianCenturies * (35999.05029 - 0.0001537 * julianCenturies);
150        }
151
152        /**
153         * Return the unitless <a href="https://en.wikipedia.org/wiki/Eccentricity_%28orbit%29">eccentricity of earth's orbit</a>.
154         * 
155         * @param julianCenturies
156         *            the number of Julian centuries since <a href=
157         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
158         * @return the unitless eccentricity
159         */
160        private static double getEarthOrbitEccentricity(double julianCenturies) {
161                return 0.016708634 - julianCenturies * (0.000042037 + 0.0000001267 * julianCenturies);
162        }
163
164        /**
165         * Returns the <a href="https://en.wikipedia.org/wiki/Equation_of_the_center">equation of center</a> for the sun in degrees.
166         * 
167         * @param julianCenturies
168         *            the number of Julian centuries since <a href=
169         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
170         * @return the equation of center for the sun in degrees
171         */
172        private static double getSunEquationOfCenter(double julianCenturies) {
173                double m = getSunGeometricMeanAnomaly(julianCenturies);
174                double mrad = Math.toRadians(m);
175                double sinm = Math.sin(mrad);
176                double sin2m = Math.sin(mrad + mrad);
177                double sin3m = Math.sin(mrad + mrad + mrad);
178                return sinm * (1.914602 - julianCenturies * (0.004817 + 0.000014 * julianCenturies)) + sin2m
179                                * (0.019993 - 0.000101 * julianCenturies) + sin3m * 0.000289;
180        }
181
182        /**
183         * Return the <a href="https://en.wikipedia.org/wiki/True_longitude">true longitude</a> of the sun.
184         * 
185         * @param julianCenturies
186         *            the number of Julian centuries since <a href=
187         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
188         * @return the sun's true longitude in degrees
189         */
190        private static double getSunTrueLongitude(double julianCenturies) {
191                double sunLongitude = getSunGeometricMeanLongitude(julianCenturies);
192                double center = getSunEquationOfCenter(julianCenturies); 
193                return sunLongitude + center;
194        }
195
196        // /**
197        // * Returns the <a href="https://en.wikipedia.org/wiki/True_anomaly">true anamoly</a> of the sun.
198        // *
199        // * @param julianCenturies
200        // * the number of Julian centuries since J2000.0
201        // * @return the sun's true anamoly in degrees
202        // */
203        // private static double getSunTrueAnomaly(double julianCenturies) {
204        // double meanAnomaly = getSunGeometricMeanAnomaly(julianCenturies);
205        // double equationOfCenter = getSunEquationOfCenter(julianCenturies);
206        //
207        // return meanAnomaly + equationOfCenter;
208        // }
209
210        /**
211         * Return the <a href="https://en.wikipedia.org/wiki/Apparent_longitude">apparent longitude</a> of the sun.
212         * 
213         * @param julianCenturies
214         *            the number of Julian centuries since <a href=
215         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
216         * @return sun's apparent longitude in degrees
217         */
218        private static double getSunApparentLongitude(double julianCenturies) {
219                double sunTrueLongitude = getSunTrueLongitude(julianCenturies);
220                double omega = 125.04 - 1934.136 * julianCenturies;
221                double lambda = sunTrueLongitude - 0.00569 - 0.00478 * Math.sin(Math.toRadians(omega));
222                return lambda;
223        }
224
225        /**
226         * Returns the mean <a href="https://en.wikipedia.org/wiki/Axial_tilt">obliquity of the ecliptic</a> (Axial tilt).
227         * 
228         * @param julianCenturies
229         *            the number of Julian centuries since <a href=
230         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
231         * @return the mean obliquity in degrees
232         */
233        private static double getMeanObliquityOfEcliptic(double julianCenturies) {
234                double seconds = 21.448 - julianCenturies
235                                * (46.8150 + julianCenturies * (0.00059 - julianCenturies * (0.001813)));
236                return 23.0 + (26.0 + (seconds / 60.0)) / 60.0;
237        }
238
239        /**
240         * Returns the corrected <a href="https://en.wikipedia.org/wiki/Axial_tilt">obliquity of the ecliptic</a> (Axial
241         * tilt).
242         * 
243         * @param julianCenturies
244         *            the number of Julian centuries since <a href=
245         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
246         * @return the corrected obliquity in degrees
247         */
248        private static double getObliquityCorrection(double julianCenturies) {
249                double obliquityOfEcliptic = getMeanObliquityOfEcliptic(julianCenturies);
250                double omega = 125.04 - 1934.136 * julianCenturies;
251                return obliquityOfEcliptic + 0.00256 * Math.cos(Math.toRadians(omega));
252        }
253
254        /**
255         * Return the <a href="https://en.wikipedia.org/wiki/Declination">declination</a> of the sun.
256         * 
257         * @param julianCenturies
258         *            the number of Julian centuries since <a href=
259         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
260         * @return
261         *            the sun's declination in degrees
262         */
263        private static double getSunDeclination(double julianCenturies) {
264                double obliquityCorrection = getObliquityCorrection(julianCenturies);
265                double lambda = getSunApparentLongitude(julianCenturies);
266                double sint = Math.sin(Math.toRadians(obliquityCorrection)) * Math.sin(Math.toRadians(lambda));
267                double theta = Math.toDegrees(Math.asin(sint));
268                return theta;
269        }
270
271        /**
272         * Return the <a href="https://en.wikipedia.org/wiki/Equation_of_time">Equation of Time</a> - the difference between
273         * true solar time and mean solar time
274         * 
275         * @param julianCenturies
276         *            the number of Julian centuries since <a href=
277         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
278         * @return equation of time in minutes of time
279         */
280        private static double getEquationOfTime(double julianCenturies) {
281                double epsilon = getObliquityCorrection(julianCenturies);
282                double geomMeanLongSun = getSunGeometricMeanLongitude(julianCenturies);
283                double eccentricityEarthOrbit = getEarthOrbitEccentricity(julianCenturies);
284                double geomMeanAnomalySun = getSunGeometricMeanAnomaly(julianCenturies);
285                double y = Math.tan(Math.toRadians(epsilon) / 2.0);
286                y *= y;
287                double sin2l0 = Math.sin(2.0 * Math.toRadians(geomMeanLongSun));
288                double sinm = Math.sin(Math.toRadians(geomMeanAnomalySun));
289                double cos2l0 = Math.cos(2.0 * Math.toRadians(geomMeanLongSun));
290                double sin4l0 = Math.sin(4.0 * Math.toRadians(geomMeanLongSun));
291                double sin2m = Math.sin(2.0 * Math.toRadians(geomMeanAnomalySun));
292                double equationOfTime = y * sin2l0 - 2.0 * eccentricityEarthOrbit * sinm + 4.0 * eccentricityEarthOrbit * y
293                                * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * eccentricityEarthOrbit * eccentricityEarthOrbit * sin2m;
294                return Math.toDegrees(equationOfTime) * 4.0;
295        }
296        
297        /**
298         * Return the <a href="https://en.wikipedia.org/wiki/Hour_angle">hour angle</a> of the sun in
299         * <a href="https://en.wikipedia.org/wiki/Radian">radians</a> at for the latitude.
300         * 
301         * @param latitude
302         *            the latitude of observer in degrees
303         * @param solarDeclination
304         *            the declination angle of sun in degrees
305         * @param zenith
306         *            the zenith
307         * @param solarEvent
308         *             If the hour angle is for {@link SolarEvent#SUNRISE SUNRISE} or {@link SolarEvent#SUNSET SUNSET}
309         * @return hour angle of sunrise in <a href="https://en.wikipedia.org/wiki/Radian">radians</a>
310         */
311        private static double getSunHourAngle(double latitude, double solarDeclination, double zenith, SolarEvent solarEvent) {
312                double latRad = Math.toRadians(latitude);
313                double sdRad = Math.toRadians(solarDeclination);
314                double hourAngle = (Math.acos(Math.cos(Math.toRadians(zenith)) / (Math.cos(latRad) * Math.cos(sdRad))
315                                - Math.tan(latRad) * Math.tan(sdRad)));
316                
317                if (solarEvent == SolarEvent.SUNSET) {
318                        hourAngle = -hourAngle;
319                }
320                return hourAngle;
321        }
322
323        /**
324         * Return the <a href="https://en.wikipedia.org/wiki/Celestial_coordinate_system">Solar Elevation</a> for the
325         * horizontal coordinate system at the given location at the given time. Can be negative if the sun is below the
326         * horizon. Not corrected for altitude.
327         * 
328         * @param calendar
329         *            time of calculation
330         * @param latitude
331         *            latitude of location for calculation
332         * @param longitude
333         *            longitude of location for calculation
334         * @return solar elevation in degrees - horizon is 0 degrees, civil twilight is -6 degrees
335         */
336
337        public static double getSolarElevation(Calendar calendar, double latitude, double longitude) {
338                double julianDay = getJulianDay(calendar);
339                double julianCenturies = getJulianCenturiesFromJulianDay(julianDay);
340                Double eot = getEquationOfTime(julianCenturies);
341                double adjustedLongitude = (calendar.get(Calendar.HOUR_OF_DAY) + 12.0)
342                                + (calendar.get(Calendar.MINUTE) + eot + calendar.get(Calendar.SECOND) / 60.0) / 60.0;
343                adjustedLongitude = -(adjustedLongitude * 360.0 / 24.0) % 360.0;
344                double hourAngle_rad = Math.toRadians(longitude - adjustedLongitude);
345                double declination = getSunDeclination(julianCenturies);
346                double dec_rad = Math.toRadians(declination);
347                double lat_rad = Math.toRadians(latitude);
348                return Math.toDegrees(Math.asin((Math.sin(lat_rad) * Math.sin(dec_rad))
349                                + (Math.cos(lat_rad) * Math.cos(dec_rad) * Math.cos(hourAngle_rad))));
350
351        }
352
353        /**
354         * Return the <a href="https://en.wikipedia.org/wiki/Celestial_coordinate_system">Solar Azimuth</a> for the
355         * horizontal coordinate system at the given location at the given time. Not corrected for altitude. True south is 0
356         * degrees.
357         * 
358         * @param calendar
359         *            time of calculation
360         * @param latitude
361         *            latitude of location for calculation
362         * @param longitude
363         *            longitude of location for calculation
364         * @return the solar azimuth
365         */
366
367        public static double getSolarAzimuth(Calendar calendar, double latitude, double longitude) {
368                double julianDay = getJulianDay(calendar);
369                double julianCenturies = getJulianCenturiesFromJulianDay(julianDay);
370                Double eot = getEquationOfTime(julianCenturies);
371                double adjustedLongitude = (calendar.get(Calendar.HOUR_OF_DAY) + 12.0)
372                                + (calendar.get(Calendar.MINUTE) + eot + calendar.get(Calendar.SECOND) / 60.0) / 60.0;
373                adjustedLongitude = -(adjustedLongitude * 360.0 / 24.0) % 360.0;
374                double hourAngle_rad = Math.toRadians(longitude - adjustedLongitude);
375                double declination = getSunDeclination(julianCenturies);
376                double dec_rad = Math.toRadians(declination);
377                double lat_rad = Math.toRadians(latitude);
378                return Math.toDegrees(Math.atan(Math.sin(hourAngle_rad)
379                                / ((Math.cos(hourAngle_rad) * Math.sin(lat_rad)) - (Math.tan(dec_rad) * Math.cos(lat_rad)))))+180;
380
381        }
382
383        /**
384         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC)
385         * of <a href="https://en.wikipedia.org/wiki/Noon#Solar_noon">solar noon</a> for the given day at the given location
386         * on earth. This implementation returns true solar noon as opposed to the time halfway between sunrise and sunset.
387         * Other calculators may return a more simplified calculation of halfway between sunrise and sunset. See <a href=
388         * "https://kosherjava.com/2020/07/02/definition-of-chatzos/">The Definition of <em>Chatzos</em></a> for details on
389         * solar noon calculations.
390         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCNoon(Calendar, GeoLocation)
391         * @see #getSolarNoonMidnightUTC(double, double, SolarEvent)
392         * 
393         * @param calendar
394         *            The Calendar representing the date to calculate solar noon for
395         * @param geoLocation
396         *            The location information used for astronomical calculating sun times. This class uses only requires
397         *            the longitude for calculating noon since it is the same time anywhere along the longitude line.
398         * @return the time in minutes from zero UTC
399         */
400        public double getUTCNoon(Calendar calendar, GeoLocation geoLocation) {
401                double noon = getSolarNoonMidnightUTC(getJulianDay(calendar), -geoLocation.getLongitude(), SolarEvent.NOON);
402                noon = noon / 60;
403                return noon > 0  ? noon % 24 : noon % 24 + 24; // ensure that the time is >= 0 and < 24
404        }
405        
406        /**
407         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a>
408         * (UTC) of the <a href="https://en.wikipedia.org/wiki/Midnight">solar midnight</a> for the end of the given civil
409         * day at the given location on earth (about 12 hours after solar noon). This implementation returns true solar
410         * midnight as opposed to the time halfway between sunrise and sunset. Other calculators may return a more
411         * simplified calculation of halfway between sunrise and sunset. See <a href=
412         * "https://kosherjava.com/2020/07/02/definition-of-chatzos/">The Definition of <em>Chatzos</em></a> for details on
413         * solar noon / midnight calculations.
414         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCNoon(Calendar, GeoLocation)
415         * @see #getSolarNoonMidnightUTC(double, double, SolarEvent)
416         * 
417         * @param calendar
418         *            The Calendar representing the date to calculate solar noon for
419         * @param geoLocation
420         *            The location information used for astronomical calculating sun times. This class uses only requires
421         *            the longitude for calculating noon since it is the same time anywhere along the longitude line.
422         * @return the time in minutes from zero UTC
423         */
424        public double getUTCMidnight(Calendar calendar, GeoLocation geoLocation) {
425                double midnight = getSolarNoonMidnightUTC(getJulianDay(calendar), -geoLocation.getLongitude(), SolarEvent.MIDNIGHT);
426                midnight = midnight / 60;
427                return midnight > 0  ? midnight % 24 : midnight % 24 + 24; // ensure that the time is >= 0 and < 24
428        }
429
430        /**
431         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC)
432         * of the current day <a href="http://en.wikipedia.org/wiki/Noon#Solar_noon">solar noon</a> or the the upcoming
433         * midnight (about 12 hours after solar noon) of the given day at the given location on earth.
434         * 
435         * @param julianDay
436         *            The Julian day since <a href=
437         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
438         * @param longitude
439         *            The longitude of observer in degrees
440         * @param solarEvent
441         *            If the calculation is for {@link SolarEvent#NOON NOON} or {@link SolarEvent#MIDNIGHT MIDNIGHT}
442         *            
443         * @return the time in minutes from zero UTC
444         * 
445         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCNoon(Calendar, GeoLocation)
446         * @see #getUTCNoon(Calendar, GeoLocation)
447         */
448        private static double getSolarNoonMidnightUTC(double julianDay, double longitude, SolarEvent solarEvent) {
449                julianDay = (solarEvent == SolarEvent.NOON) ? julianDay : julianDay + 0.5;
450                // First pass for approximate solar noon to calculate equation of time
451                double tnoon = getJulianCenturiesFromJulianDay(julianDay + longitude / 360.0);
452                double equationOfTime = getEquationOfTime(tnoon);
453                double solNoonUTC = (longitude * 4) - equationOfTime; // minutes
454                
455                // second pass
456                double newt = getJulianCenturiesFromJulianDay(julianDay + solNoonUTC / 1440.0);
457                equationOfTime = getEquationOfTime(newt);
458                return (solarEvent == SolarEvent.NOON ? 720 : 1440) + (longitude * 4) - equationOfTime;
459        }
460        
461        /**
462         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC)
463         * of sunrise or sunset in minutes for the given day at the given location on earth.
464         * @todo Possibly increase the number of passes for improved accuracy, especially in the Arctic areas.
465         * 
466         * @param calendar
467         *            The calendar
468         * @param latitude
469         *            The latitude of observer in degrees
470         * @param longitude
471         *            Longitude of observer in degrees
472         * @param zenith
473         *            Zenith
474         * @param solarEvent
475         *             If the calculation is for {@link SolarEvent#SUNRISE SUNRISE} or {@link SolarEvent#SUNSET SUNSET}
476         * @return the time in minutes from zero Universal Coordinated Time (UTC)
477         */
478        private static double getSunRiseSetUTC(Calendar calendar, double latitude, double longitude, double zenith,
479                        SolarEvent solarEvent) {
480                double julianDay = getJulianDay(calendar);
481
482                // Find the time of solar noon at the location, and use that declination.
483                // This is better than start of the Julian day
484                // TODO really not needed since the Julian day starts from local fixed noon. Changing this would be more
485                // efficient but would likely cause a very minor discrepancy in the calculated times (likely not reducing
486                // accuracy, just slightly different, thus potentially breaking test cases). Regardless, it would be within
487                // milliseconds.
488                double noonmin = getSolarNoonMidnightUTC(julianDay, longitude, SolarEvent.NOON);
489                                                                                                                                                                                
490                double tnoon = getJulianCenturiesFromJulianDay(julianDay + noonmin / 1440.0);
491
492                // First calculates sunrise and approximate length of day
493                double equationOfTime = getEquationOfTime(tnoon);
494                double solarDeclination = getSunDeclination(tnoon);
495                double hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
496                double delta = longitude - Math.toDegrees(hourAngle);
497                double timeDiff = 4 * delta;
498                double timeUTC = 720 + timeDiff - equationOfTime;
499
500                // Second pass includes fractional Julian Day in gamma calc
501                double newt = getJulianCenturiesFromJulianDay(julianDay + timeUTC / 1440.0);
502                equationOfTime = getEquationOfTime(newt);
503                
504                solarDeclination = getSunDeclination(newt);
505                hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
506                delta = longitude - Math.toDegrees(hourAngle);
507                timeDiff = 4 * delta;
508                timeUTC = 720 + timeDiff - equationOfTime;
509                return timeUTC;
510        }
511}