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 the start of the Julian
100         *         day. Fractional days / time 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         * Return the <a href="https://en.wikipedia.org/wiki/Apparent_longitude">apparent longitude</a> of the sun.
198         * 
199         * @param julianCenturies
200         *            the number of Julian centuries since <a href=
201         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
202         * @return sun's apparent longitude in degrees
203         */
204        private static double getSunApparentLongitude(double julianCenturies) {
205                double sunTrueLongitude = getSunTrueLongitude(julianCenturies);
206                double omega = 125.04 - 1934.136 * julianCenturies;
207                double lambda = sunTrueLongitude - 0.00569 - 0.00478 * Math.sin(Math.toRadians(omega));
208                return lambda;
209        }
210
211        /**
212         * Returns the mean <a href="https://en.wikipedia.org/wiki/Axial_tilt">obliquity of the ecliptic</a> (Axial tilt).
213         * 
214         * @param julianCenturies
215         *            the number of Julian centuries since <a href=
216         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
217         * @return the mean obliquity in degrees
218         */
219        private static double getMeanObliquityOfEcliptic(double julianCenturies) {
220                double seconds = 21.448 - julianCenturies
221                                * (46.8150 + julianCenturies * (0.00059 - julianCenturies * (0.001813)));
222                return 23.0 + (26.0 + (seconds / 60.0)) / 60.0;
223        }
224
225        /**
226         * Returns the corrected <a href="https://en.wikipedia.org/wiki/Axial_tilt">obliquity of the ecliptic</a> (Axial
227         * tilt).
228         * 
229         * @param julianCenturies
230         *            the number of Julian centuries since <a href=
231         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
232         * @return the corrected obliquity in degrees
233         */
234        private static double getObliquityCorrection(double julianCenturies) {
235                double obliquityOfEcliptic = getMeanObliquityOfEcliptic(julianCenturies);
236                double omega = 125.04 - 1934.136 * julianCenturies;
237                return obliquityOfEcliptic + 0.00256 * Math.cos(Math.toRadians(omega));
238        }
239
240        /**
241         * Return the <a href="https://en.wikipedia.org/wiki/Declination">declination</a> of the sun.
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
247         *            the sun's declination in degrees
248         */
249        private static double getSunDeclination(double julianCenturies) {
250                double obliquityCorrection = getObliquityCorrection(julianCenturies);
251                double lambda = getSunApparentLongitude(julianCenturies);
252                double sint = Math.sin(Math.toRadians(obliquityCorrection)) * Math.sin(Math.toRadians(lambda));
253                double theta = Math.toDegrees(Math.asin(sint));
254                return theta;
255        }
256
257        /**
258         * Return the <a href="https://en.wikipedia.org/wiki/Equation_of_time">Equation of Time</a> - the difference between
259         * true solar time and mean solar time
260         * 
261         * @param julianCenturies
262         *            the number of Julian centuries since <a href=
263         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
264         * @return equation of time in minutes of time
265         */
266        private static double getEquationOfTime(double julianCenturies) {
267                double epsilon = getObliquityCorrection(julianCenturies);
268                double geomMeanLongSun = getSunGeometricMeanLongitude(julianCenturies);
269                double eccentricityEarthOrbit = getEarthOrbitEccentricity(julianCenturies);
270                double geomMeanAnomalySun = getSunGeometricMeanAnomaly(julianCenturies);
271                double y = Math.tan(Math.toRadians(epsilon) / 2.0);
272                y *= y;
273                double sin2l0 = Math.sin(2.0 * Math.toRadians(geomMeanLongSun));
274                double sinm = Math.sin(Math.toRadians(geomMeanAnomalySun));
275                double cos2l0 = Math.cos(2.0 * Math.toRadians(geomMeanLongSun));
276                double sin4l0 = Math.sin(4.0 * Math.toRadians(geomMeanLongSun));
277                double sin2m = Math.sin(2.0 * Math.toRadians(geomMeanAnomalySun));
278                double equationOfTime = y * sin2l0 - 2.0 * eccentricityEarthOrbit * sinm + 4.0 * eccentricityEarthOrbit * y
279                                * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * eccentricityEarthOrbit * eccentricityEarthOrbit * sin2m;
280                return Math.toDegrees(equationOfTime) * 4.0;
281        }
282        
283        /**
284         * Return the <a href="https://en.wikipedia.org/wiki/Hour_angle">hour angle</a> of the sun in
285         * <a href="https://en.wikipedia.org/wiki/Radian">radians</a> at for the latitude.
286         * 
287         * @param latitude
288         *            the latitude of observer in degrees
289         * @param solarDeclination
290         *            the declination angle of sun in degrees
291         * @param zenith
292         *            the zenith
293         * @param solarEvent
294         *             If the hour angle is for {@link SolarEvent#SUNRISE SUNRISE} or {@link SolarEvent#SUNSET SUNSET}
295         * @return hour angle of sunrise in <a href="https://en.wikipedia.org/wiki/Radian">radians</a>
296         */
297        private static double getSunHourAngle(double latitude, double solarDeclination, double zenith, SolarEvent solarEvent) {
298                double latRad = Math.toRadians(latitude);
299                double sdRad = Math.toRadians(solarDeclination);
300                double hourAngle = (Math.acos(Math.cos(Math.toRadians(zenith)) / (Math.cos(latRad) * Math.cos(sdRad))
301                                - Math.tan(latRad) * Math.tan(sdRad)));
302                
303                if (solarEvent == SolarEvent.SUNSET) {
304                        hourAngle = -hourAngle;
305                }
306                return hourAngle;
307        }
308        
309        /**
310         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getSolarElevation(Calendar, GeoLocation)
311         */
312        public double getSolarElevation(Calendar calendar, GeoLocation geoLocation) {
313                return getSolarElevationAzimuth(calendar, geoLocation, false);
314
315        }
316        
317        /**
318         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getSolarAzimuth(Calendar, GeoLocation)
319         */
320        public double getSolarAzimuth(Calendar calendar, GeoLocation geoLocation) {
321                return getSolarElevationAzimuth(calendar, geoLocation, true);
322        }
323
324        /**
325         * Return the <a href="https://en.wikipedia.org/wiki/Celestial_coordinate_system">Solar Elevation</a> or
326         * <a href="https://en.wikipedia.org/wiki/Celestial_coordinate_system">Solar Azimuth</a> at the given location
327         * and time. Can be negative if the sun is below the horizon. Elevation is based on sea-level and is not
328         * adjusted for altitude.
329         * 
330         * @param calendar
331         *            time of calculation
332         * @param geoLocation
333         *            The location for calculating the elevation or azimuth.
334         * @param isAzimuth
335         *            true for azimuth, false for elevation
336         * @return solar elevation or azimuth in degrees.
337         * 
338         * @see #getSolarElevation(Calendar, GeoLocation)
339         * @see #getSolarAzimuth(Calendar, GeoLocation)
340         */
341        private double getSolarElevationAzimuth(Calendar calendar, GeoLocation geoLocation, boolean isAzimuth) {
342                double latitude = geoLocation.getLatitude();
343                double longitude = geoLocation.getLongitude();
344                
345                Calendar cloned = (Calendar) calendar.clone();
346                int offset = - adjustHourForTimeZone(cloned);
347                cloned.add(Calendar.MILLISECOND, offset);
348                int minute = cloned.get(Calendar.MINUTE);
349                int second = cloned.get(Calendar.SECOND);
350                int hour = cloned.get(Calendar.HOUR_OF_DAY);
351                int milli = cloned.get(Calendar.MILLISECOND);
352                
353                double time = (hour + (minute + (second + (milli / 1000.0)) / 60.0) / 60.0 ) / 24.0;
354                double julianDay = getJulianDay(cloned) + time;
355                double julianCenturies = getJulianCenturiesFromJulianDay(julianDay);
356                double eot = getEquationOfTime(julianCenturies);
357                double theta = getSunDeclination(julianCenturies);
358                
359                double adjustment = time + eot / 1440;
360                double trueSolarTime = ((adjustment + longitude / 360) + 2) % 1; // adding 2 to ensure that it never ends up negative
361                double hourAngelRad = trueSolarTime * Math.PI * 2 - Math.PI;
362                double cosZenith = Math.sin(Math.toRadians(latitude)) * Math.sin(Math.toRadians(theta))
363                                +  Math.cos(Math.toRadians(latitude)) * Math.cos(Math.toRadians(theta)) * Math.cos(hourAngelRad);
364                double zenith = Math.toDegrees(Math.acos(cosZenith > 1 ? 1 : cosZenith < -1 ? -1 : cosZenith));
365                double azDenom = Math.cos(Math.toRadians(latitude)) * Math.sin(Math.toRadians(zenith));
366                double refractionAdjustment = 0;
367                double elevation = 90.0 - (zenith - refractionAdjustment);
368                double azimuth = 0;
369                double azRad = (Math.sin(Math.toRadians(latitude)) * Math.cos(Math.toRadians(zenith))
370                                - Math.sin(Math.toRadians(theta))) / azDenom;
371                if(Math.abs(azDenom) > 0.001) {
372                        azimuth = 180 - Math.toDegrees(Math.acos(azRad > 1 ? 1 : azRad < -1? -1 : azRad)) * (hourAngelRad > 0 ? -1 : 1) ;
373                } else {
374                        azimuth = latitude > 0 ? 180 : 0;
375                }
376                return isAzimuth ? azimuth % 360 : elevation;
377        }
378        
379        /**
380         * Returns the hour of day adjusted for the timezone and DST. This is needed for the azimuth and elevation
381         * calculations.
382         * @param calendar the Calendar to extract the hour from. This must have the timezone set to the proper timezone.
383         * @return the adjusted hour corrected for timezone and DST offset.
384         */
385        private int adjustHourForTimeZone(Calendar calendar) {
386                int offset = calendar.getTimeZone().getRawOffset();
387                int dstOffset = calendar.getTimeZone().getDSTSavings();
388                if(calendar.getTimeZone().inDaylightTime(calendar.getTime())) {
389                        offset = offset + dstOffset;
390                }
391                return offset;
392        }
393
394        /**
395         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC)
396         * of <a href="https://en.wikipedia.org/wiki/Noon#Solar_noon">solar noon</a> for the given day at the given location
397         * on earth. This implementation returns true solar noon as opposed to the time halfway between sunrise and sunset.
398         * Other calculators may return a more simplified calculation of halfway between sunrise and sunset. See <a href=
399         * "https://kosherjava.com/2020/07/02/definition-of-chatzos/">The Definition of <em>Chatzos</em></a> for details on
400         * solar noon calculations.
401         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCNoon(Calendar, GeoLocation)
402         * @see #getSolarNoonMidnightUTC(double, double, SolarEvent)
403         * 
404         * @param calendar
405         *            The Calendar representing the date to calculate solar noon for
406         * @param geoLocation
407         *            The location information used for astronomical calculating sun times. This class uses only requires
408         *            the longitude for calculating noon since it is the same time anywhere along the longitude line.
409         * @return the time in minutes from zero UTC
410         */
411        public double getUTCNoon(Calendar calendar, GeoLocation geoLocation) {
412                double noon = getSolarNoonMidnightUTC(getJulianDay(calendar), -geoLocation.getLongitude(), SolarEvent.NOON);
413                noon = noon / 60;
414                return noon > 0  ? noon % 24 : noon % 24 + 24; // ensure that the time is >= 0 and < 24
415        }
416        
417        /**
418         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a>
419         * (UTC) of the <a href="https://en.wikipedia.org/wiki/Midnight">solar midnight</a> for the end of the given civil
420         * day at the given location on earth (about 12 hours after solar noon). This implementation returns true solar
421         * midnight as opposed to the time halfway between sunrise and sunset. Other calculators may return a more
422         * simplified calculation of halfway between sunrise and sunset. See <a href=
423         * "https://kosherjava.com/2020/07/02/definition-of-chatzos/">The Definition of <em>Chatzos</em></a> for details on
424         * solar noon / midnight calculations.
425         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCNoon(Calendar, GeoLocation)
426         * @see #getSolarNoonMidnightUTC(double, double, SolarEvent)
427         * 
428         * @param calendar
429         *            The Calendar representing the date to calculate solar noon for
430         * @param geoLocation
431         *            The location information used for astronomical calculating sun times. This class uses only requires
432         *            the longitude for calculating noon since it is the same time anywhere along the longitude line.
433         * @return the time in minutes from zero UTC
434         */
435        public double getUTCMidnight(Calendar calendar, GeoLocation geoLocation) {
436                double midnight = getSolarNoonMidnightUTC(getJulianDay(calendar), -geoLocation.getLongitude(), SolarEvent.MIDNIGHT);
437                midnight = midnight / 60;
438                return midnight > 0  ? midnight % 24 : midnight % 24 + 24; // ensure that the time is >= 0 and < 24
439        }
440
441        /**
442         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC)
443         * of the current day <a href="http://en.wikipedia.org/wiki/Noon#Solar_noon">solar noon</a> or the the upcoming
444         * midnight (about 12 hours after solar noon) of the given day at the given location on earth.
445         * 
446         * @param julianDay
447         *            The Julian day since <a href=
448         *            "https://en.wikipedia.org/wiki/Epoch_(astronomy)#J2000">J2000.0</a>.
449         * @param longitude
450         *            The longitude of observer in degrees
451         * @param solarEvent
452         *            If the calculation is for {@link SolarEvent#NOON NOON} or {@link SolarEvent#MIDNIGHT MIDNIGHT}
453         *            
454         * @return the time in minutes from zero UTC
455         * 
456         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCNoon(Calendar, GeoLocation)
457         * @see #getUTCNoon(Calendar, GeoLocation)
458         */
459        private static double getSolarNoonMidnightUTC(double julianDay, double longitude, SolarEvent solarEvent) {
460                julianDay = (solarEvent == SolarEvent.NOON) ? julianDay : julianDay + 0.5;
461                // First pass for approximate solar noon to calculate equation of time
462                double tnoon = getJulianCenturiesFromJulianDay(julianDay + longitude / 360.0);
463                double equationOfTime = getEquationOfTime(tnoon);
464                double solNoonUTC = (longitude * 4) - equationOfTime; // minutes
465                
466                // second pass
467                double newt = getJulianCenturiesFromJulianDay(julianDay + solNoonUTC / 1440.0);
468                equationOfTime = getEquationOfTime(newt);
469                return (solarEvent == SolarEvent.NOON ? 720 : 1440) + (longitude * 4) - equationOfTime;
470        }
471        
472        /**
473         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC)
474         * of sunrise or sunset in minutes for the given day at the given location on earth.
475         * @todo Possibly increase the number of passes for improved accuracy, especially in the Arctic areas.
476         * 
477         * @param calendar
478         *            The calendar
479         * @param latitude
480         *            The latitude of observer in degrees
481         * @param longitude
482         *            Longitude of observer in degrees
483         * @param zenith
484         *            Zenith
485         * @param solarEvent
486         *             If the calculation is for {@link SolarEvent#SUNRISE SUNRISE} or {@link SolarEvent#SUNSET SUNSET}
487         * @return the time in minutes from zero Universal Coordinated Time (UTC)
488         */
489        private static double getSunRiseSetUTC(Calendar calendar, double latitude, double longitude, double zenith,
490                        SolarEvent solarEvent) {
491                double julianDay = getJulianDay(calendar);
492
493                // Find the time of solar noon at the location, and use that declination.
494                // This is better than start of the Julian day
495                // TODO really not needed since the Julian day starts from local fixed noon. Changing this would be more
496                // efficient but would likely cause a very minor discrepancy in the calculated times (likely not reducing
497                // accuracy, just slightly different, thus potentially breaking test cases). Regardless, it would be within
498                // milliseconds.
499                double noonmin = getSolarNoonMidnightUTC(julianDay, longitude, SolarEvent.NOON);
500                                                                                                                                                                                
501                double tnoon = getJulianCenturiesFromJulianDay(julianDay + noonmin / 1440.0);
502
503                // First calculates sunrise and approximate length of day
504                double equationOfTime = getEquationOfTime(tnoon);
505                double solarDeclination = getSunDeclination(tnoon);
506                double hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
507                double delta = longitude - Math.toDegrees(hourAngle);
508                double timeDiff = 4 * delta;
509                double timeUTC = 720 + timeDiff - equationOfTime;
510
511                // Second pass includes fractional Julian Day in gamma calc
512                double newt = getJulianCenturiesFromJulianDay(julianDay + timeUTC / 1440.0);
513                equationOfTime = getEquationOfTime(newt);
514                
515                solarDeclination = getSunDeclination(newt);
516                hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
517                delta = longitude - Math.toDegrees(hourAngle);
518                timeDiff = 4 * delta;
519                timeUTC = 720 + timeDiff - equationOfTime;
520                return timeUTC;
521        }
522}