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. This calculator uses the Java algorithm
022 * written by <a href="https://web.archive.org/web/20090531215353/http://www.kevinboone.com/suntimes.html">Kevin
023 * Boone</a> that is based on the <a href = "https://aa.usno.navy.mil/">US Naval Observatory's</a><a
024 * href="https://aa.usno.navy.mil/publications/asa">Astronomical Almanac</a> and used with his permission. Added to Kevin's
025 * code is adjustment of the zenith to account for elevation. This algorithm returns the same time every year and does not
026 * account for leap years. It is not as accurate as the Jean Meeus based {@link NOAACalculator} that is the default calculator
027 * use by the KosherJava <em>zmanim</em> library.
028 *
029 * @author &copy; Eliyahu Hershfeld 2004 - 2025
030 * @author &copy; Kevin Boone 2000
031 */
032public class SunTimesCalculator extends AstronomicalCalculator {
033        
034        /**
035         * Default constructor of the SunTimesCalculator.
036         */
037        public SunTimesCalculator() {
038                super();
039        }
040
041        /**
042         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getCalculatorName()
043         */
044        public String getCalculatorName() {
045                return "US Naval Almanac Algorithm";
046        }
047
048        /**
049         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCSunrise(Calendar, GeoLocation, double, boolean)
050         */
051        public double getUTCSunrise(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) {
052                double elevation = adjustForElevation ? geoLocation.getElevation() : 0;
053                double adjustedZenith = adjustZenith(zenith, elevation);
054                return getTimeUTC(calendar, geoLocation, adjustedZenith, true);
055        }
056
057        /**
058         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCSunset(Calendar, GeoLocation, double, boolean)
059         */
060        public double getUTCSunset(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) {
061                double elevation = adjustForElevation ? geoLocation.getElevation() : 0;
062                double adjustedZenith = adjustZenith(zenith, elevation);
063                return getTimeUTC(calendar, geoLocation, adjustedZenith, false);
064        }
065
066        /**
067         * The number of degrees of longitude that corresponds to one hour of time difference.
068         */
069        private static final double DEG_PER_HOUR = 360.0 / 24.0;
070
071        /**
072         * The sine in degrees.
073         * @param deg the degrees
074         * @return sin of the angle in degrees
075         */
076        private static double sinDeg(double deg) {
077                return Math.sin(deg * 2.0 * Math.PI / 360.0);
078        }
079
080        /**
081         * Return the arc cosine in degrees.
082         * @param x angle
083         * @return acos of the angle in degrees
084         */
085        private static double acosDeg(double x) {
086                return Math.acos(x) * 360.0 / (2 * Math.PI);
087        }
088
089        /**
090         * Return the arc sine in degrees.
091         * @param x angle
092         * @return asin of the angle in degrees
093         */
094        private static double asinDeg(double x) {
095                return Math.asin(x) * 360.0 / (2 * Math.PI);
096        }
097
098        /**
099         * Return the tangent in degrees.
100         * @param deg degrees
101         * @return tan of the angle in degrees
102         */
103        private static double tanDeg(double deg) {
104                return Math.tan(deg * 2.0 * Math.PI / 360.0);
105        }
106        
107        /**
108         * Calculate cosine of the angle in degrees
109         * 
110         * @param deg degrees
111         * @return cosine of the angle in degrees
112         */
113        private static double cosDeg(double deg) {
114                return Math.cos(deg * 2.0 * Math.PI / 360.0);
115        }
116
117        /**
118         * Get time difference between location's longitude and the Meridian, in hours.
119         * 
120         * @param longitude the longitude
121         * @return time difference between the location's longitude and the Meridian, in hours. West of Meridian has a negative time difference
122         */
123        private static double getHoursFromMeridian(double longitude) {
124                return longitude / DEG_PER_HOUR;
125        }
126        
127        /**
128         * Calculate the approximate time of sunset or sunrise in days since midnight Jan 1st, assuming 6am and 6pm events. We
129         * need this figure to derive the Sun's mean anomaly.
130         * 
131         * @param dayOfYear the day of year
132         * @param hoursFromMeridian hours from the meridian
133         * @param isSunrise true for sunrise and false for sunset
134         * 
135         * @return the approximate time of sunset or sunrise in days since midnight Jan 1st, assuming 6am and 6pm events. We
136         * need this figure to derive the Sun's mean anomaly.
137         */
138        private static double getApproxTimeDays(int dayOfYear, double hoursFromMeridian, boolean isSunrise) {
139                if (isSunrise) {
140                        return dayOfYear + ((6.0 - hoursFromMeridian) / 24);
141                } else { // sunset
142                        return dayOfYear + ((18.0 - hoursFromMeridian) / 24);
143                }
144        }
145        
146        /**
147         * Calculate the Sun's mean anomaly in degrees, at sunrise or sunset, given the longitude in degrees
148         * 
149         * @param dayOfYear the day of the year
150         * @param longitude longitude
151         * @param isSunrise true for sunrise and false for sunset
152         * @return the Sun's mean anomaly in degrees
153         */
154        private static double getMeanAnomaly(int dayOfYear, double longitude, boolean isSunrise) {
155                return (0.9856 * getApproxTimeDays(dayOfYear, getHoursFromMeridian(longitude), isSunrise)) - 3.289;
156        }
157
158        /**
159         * Returns the Sun's true longitude in degrees. 
160         * @param sunMeanAnomaly the Sun's mean anomaly in degrees
161         * @return the Sun's true longitude in degrees. The result is an angle &gt;= 0 and &lt;= 360.
162         */
163        private static double getSunTrueLongitude(double sunMeanAnomaly) {
164                double l = sunMeanAnomaly + (1.916 * sinDeg(sunMeanAnomaly)) + (0.020 * sinDeg(2 * sunMeanAnomaly)) + 282.634;
165
166                // get longitude into 0-360 degree range
167                if (l >= 360.0) {
168                        l = l - 360.0;
169                }
170                if (l < 0) {
171                        l = l + 360.0;
172                }
173                return l;
174        }
175
176        /**
177         * Calculates the Sun's right ascension in hours.
178         * @param sunTrueLongitude the Sun's true longitude in degrees &gt; 0 and &lt; 360.
179         * @return the Sun's right ascension in hours in angles &gt; 0 and &lt; 360.
180         */
181        private static double getSunRightAscensionHours(double sunTrueLongitude) {
182                double a = 0.91764 * tanDeg(sunTrueLongitude);
183                double ra = 360.0 / (2.0 * Math.PI) * Math.atan(a);
184
185                double lQuadrant = Math.floor(sunTrueLongitude / 90.0) * 90.0;
186                double raQuadrant = Math.floor(ra / 90.0) * 90.0;
187                ra = ra + (lQuadrant - raQuadrant);
188
189                return ra / DEG_PER_HOUR; // convert to hours
190        }
191
192        /**
193         * Calculate the cosine of the Sun's local hour angle
194         * 
195         * @param sunTrueLongitude the sun's true longitude
196         * @param latitude the latitude
197         * @param zenith the zenith
198         * @return the cosine of the Sun's local hour angle
199         */
200        private static double getCosLocalHourAngle(double sunTrueLongitude, double latitude, double zenith) {
201                double sinDec = 0.39782 * sinDeg(sunTrueLongitude);
202                double cosDec = cosDeg(asinDeg(sinDec));
203                return (cosDeg(zenith) - (sinDec * sinDeg(latitude))) / (cosDec * cosDeg(latitude));
204        }
205        
206        /**
207         * Calculate local mean time of rising or setting. By 'local' is meant the exact time at the location, assuming that
208         * there were no time zone. That is, the time difference between the location and the Meridian depended entirely on
209         * the longitude. We can't do anything with this time directly; we must convert it to UTC and then to a local time.
210         * 
211         * @param localHour the local hour
212         * @param sunRightAscensionHours the sun's right ascension in hours
213         * @param approxTimeDays approximate time days
214         * 
215         * @return the fractional number of hours since midnight as a double
216         */
217        private static double getLocalMeanTime(double localHour, double sunRightAscensionHours, double approxTimeDays) {
218                return localHour + sunRightAscensionHours - (0.06571 * approxTimeDays) - 6.622;
219        }
220
221        /**
222         * Get sunrise or sunset time in UTC, according to flag. This time is returned as
223         * a double and is not adjusted for time-zone.
224         * 
225         * @param calendar
226         *            the Calendar object to extract the day of year for calculation
227         * @param geoLocation
228         *            the GeoLocation object that contains the latitude and longitude
229         * @param zenith
230         *            Sun's zenith, in degrees
231         * @param isSunrise
232         *            True for sunrise and false for sunset.
233         * @return the time as a double. If an error was encountered in the calculation
234         *         (expected behavior for some locations such as near the poles,
235         *         {@link Double#NaN} will be returned.
236         */
237        private static double getTimeUTC(Calendar calendar, GeoLocation geoLocation, double zenith, boolean isSunrise) {
238                int dayOfYear = calendar.get(Calendar.DAY_OF_YEAR);
239                double sunMeanAnomaly = getMeanAnomaly(dayOfYear, geoLocation.getLongitude(), isSunrise);
240                double sunTrueLong = getSunTrueLongitude(sunMeanAnomaly);
241                double sunRightAscensionHours = getSunRightAscensionHours(sunTrueLong);
242                double cosLocalHourAngle = getCosLocalHourAngle(sunTrueLong, geoLocation.getLatitude(), zenith);
243
244                double localHourAngle;
245                if (isSunrise) {
246                        localHourAngle = 360.0 - acosDeg(cosLocalHourAngle);
247                } else { // sunset
248                        localHourAngle = acosDeg(cosLocalHourAngle);
249                }
250                double localHour = localHourAngle / DEG_PER_HOUR;
251
252                double localMeanTime = getLocalMeanTime(localHour, sunRightAscensionHours,
253                                getApproxTimeDays(dayOfYear, getHoursFromMeridian(geoLocation.getLongitude()), isSunrise));
254                double pocessedTime = localMeanTime - getHoursFromMeridian(geoLocation.getLongitude());
255                return pocessedTime > 0  ? pocessedTime % 24 : pocessedTime % 24 + 24; // ensure that the time is >= 0 and < 24
256        }
257        
258        /**
259         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC)
260         * of <a href="https://en.wikipedia.org/wiki/Noon#Solar_noon">solar noon</a> for the given day at the given location
261         * on earth. This implementation returns solar noon as the time halfway between sunrise and sunset.
262         * {@link NOAACalculator}, the default calculator, returns true solar noon. See <a href=
263         * "https://kosherjava.com/2020/07/02/definition-of-chatzos/">The Definition of Chatzos</a> for details on solar
264         * noon calculations.
265         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCNoon(Calendar, GeoLocation)
266         * @see NOAACalculator
267         * 
268         * @param calendar
269         *            The Calendar representing the date to calculate solar noon for
270         * @param geoLocation
271         *            The location information used for astronomical calculating sun times.
272         * @return the time in minutes from zero UTC. If an error was encountered in the calculation (expected behavior for
273         *         some locations such as near the poles, {@link Double#NaN} will be returned.
274         */
275        public double getUTCNoon(Calendar calendar, GeoLocation geoLocation) {
276                double sunrise = getUTCSunrise(calendar, geoLocation, 90, false);
277                double sunset = getUTCSunset(calendar, geoLocation, 90, false);
278                double noon = sunrise + ((sunset - sunrise) / 2);
279                if (noon < 0) {
280                        noon += 12;
281                }
282                if (noon < sunrise) {
283                        noon -= 12;
284                } 
285                return noon;
286        }
287        
288        /**
289         * Return the <a href="https://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC)
290         * of midnight for the given day at the given location on earth. This implementation returns solar midnight as 12 hours
291         * after utc noon that is  halfway between sunrise and sunset.
292         * {@link NOAACalculator}, the default calculator, returns true solar noon. See <a href=
293         * "https://kosherjava.com/2020/07/02/definition-of-chatzos/">The Definition of Chatzos</a> for details on solar
294         * noon calculations.
295         * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getUTCNoon(Calendar, GeoLocation)
296         * @see NOAACalculator
297         * 
298         * @param calendar
299         *            The Calendar representing the date to calculate solar noon for
300         * @param geoLocation
301         *            The location information used for astronomical calculating sun times.
302         * @return the time in minutes from zero UTC. If an error was encountered in the calculation (expected behavior for
303         *         some locations such as near the poles, {@link Double#NaN} will be returned.
304         */
305        public double getUTCMidnight(Calendar calendar, GeoLocation geoLocation) {
306                return (getUTCNoon(calendar, geoLocation) + 12);
307        }
308}