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