diff --git a/src/main/java/com/kosherjava/zmanim/util/NOAACalculator.java b/src/main/java/com/kosherjava/zmanim/util/NOAACalculator.java index 799c3a6e..f615bd2c 100644 --- a/src/main/java/com/kosherjava/zmanim/util/NOAACalculator.java +++ b/src/main/java/com/kosherjava/zmanim/util/NOAACalculator.java @@ -1,6 +1,6 @@ /* * Zmanim Java API - * Copyright (C) 2004-2023 Eliyahu Hershfeld + * Copyright (C) 2004-2024 Eliyahu Hershfeld * * This library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General * Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) @@ -17,6 +17,7 @@ import java.util.Calendar; + /** * Implementation of sunrise and sunset methods to calculate astronomical times based on the NOAA algorithm. This calculator uses the Java algorithm based on the implementation by Wikipedia Sunrise Equation article. * - * @author © Eliyahu Hershfeld 2011 - 2023 + * @author © Eliyahu Hershfeld 2011 - 2024 */ public class NOAACalculator extends AstronomicalCalculator { /** @@ -41,12 +42,21 @@ public class NOAACalculator extends AstronomicalCalculator { * Julian days per century. */ private static final double JULIAN_DAYS_PER_CENTURY = 36525.0; + + /** + * An enum to indicate what type of solar event is being calculated. + */ + protected enum SolarEvent { + /**SUNRISE A solar event related to sunrise*/SUNRISE, /**SUNSET A solar event related to sunset*/SUNSET + // possibly add the following in the future, if added, an IllegalArgumentException should be thrown in getSunHourAngle + // /**NOON A solar event related to noon*/NOON, /**MIDNIGHT A solar event related to midnight*/MIDNIGHT + } /** * @see com.kosherjava.zmanim.util.AstronomicalCalculator#getCalculatorName() */ public String getCalculatorName() { - return "US National Oceanic and Atmospheric Administration Algorithm"; + return "US National Oceanic and Atmospheric Administration Algorithm"; // Implementation of the Jean Meeus algorithm } /** @@ -55,19 +65,11 @@ public String getCalculatorName() { public double getUTCSunrise(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) { double elevation = adjustForElevation ? geoLocation.getElevation() : 0; double adjustedZenith = adjustZenith(zenith, elevation); - double sunrise = getSunriseUTC(getJulianDay(calendar), geoLocation.getLatitude(), -geoLocation.getLongitude(), adjustedZenith); sunrise = sunrise / 60; - // ensure that the time is >= 0 and < 24 - while (sunrise < 0.0) { - sunrise += 24.0; - } - while (sunrise >= 24.0) { - sunrise -= 24.0; - } - return sunrise; + return sunrise > 0 ? sunrise % 24 : sunrise % 24 + 24; // ensure that the time is >= 0 and < 24 } /** @@ -76,19 +78,10 @@ public double getUTCSunrise(Calendar calendar, GeoLocation geoLocation, double z public double getUTCSunset(Calendar calendar, GeoLocation geoLocation, double zenith, boolean adjustForElevation) { double elevation = adjustForElevation ? geoLocation.getElevation() : 0; double adjustedZenith = adjustZenith(zenith, elevation); - double sunset = getSunsetUTC(getJulianDay(calendar), geoLocation.getLatitude(), -geoLocation.getLongitude(), adjustedZenith); sunset = sunset / 60; - - // ensure that the time is >= 0 and < 24 - while (sunset < 0.0) { - sunset += 24.0; - } - while (sunset >= 24.0) { - sunset -= 24.0; - } - return sunset; + return sunset > 0 ? sunset % 24 : sunset % 24 + 24; // ensure that the time is >= 0 and < 24 } /** @@ -109,7 +102,6 @@ private static double getJulianDay(Calendar calendar) { } int a = year / 100; int b = 2 - a + a / 4; - return Math.floor(365.25 * (year + 4716)) + Math.floor(30.6001 * (month + 1)) + day + b - 1524.5; } @@ -148,14 +140,7 @@ private static double getJulianDayFromJulianCenturies(double julianCenturies) { */ private static double getSunGeometricMeanLongitude(double julianCenturies) { double longitude = 280.46646 + julianCenturies * (36000.76983 + 0.0003032 * julianCenturies); - while (longitude > 360.0) { - longitude -= 360.0; - } - while (longitude < 0.0) { - longitude += 360.0; - } - - return longitude; // in degrees + return longitude > 0 ? longitude % 360 : longitude % 360 + 360; // ensure that the longitude is >= 0 and < 360 } /** @@ -192,14 +177,12 @@ private static double getEarthOrbitEccentricity(double julianCenturies) { */ private static double getSunEquationOfCenter(double julianCenturies) { double m = getSunGeometricMeanAnomaly(julianCenturies); - double mrad = Math.toRadians(m); double sinm = Math.sin(mrad); double sin2m = Math.sin(mrad + mrad); double sin3m = Math.sin(mrad + mrad + mrad); - return sinm * (1.914602 - julianCenturies * (0.004817 + 0.000014 * julianCenturies)) + sin2m - * (0.019993 - 0.000101 * julianCenturies) + sin3m * 0.000289;// in degrees + * (0.019993 - 0.000101 * julianCenturies) + sin3m * 0.000289; // in degrees } /** @@ -212,8 +195,7 @@ private static double getSunEquationOfCenter(double julianCenturies) { */ private static double getSunTrueLongitude(double julianCenturies) { double sunLongitude = getSunGeometricMeanLongitude(julianCenturies); - double center = getSunEquationOfCenter(julianCenturies); - + double center = getSunEquationOfCenter(julianCenturies); return sunLongitude + center; // in degrees } @@ -241,9 +223,9 @@ private static double getSunTrueLongitude(double julianCenturies) { */ private static double getSunApparentLongitude(double julianCenturies) { double sunTrueLongitude = getSunTrueLongitude(julianCenturies); - double omega = 125.04 - 1934.136 * julianCenturies; - return sunTrueLongitude - 0.00569 - 0.00478 * Math.sin(Math.toRadians(omega)); // in degrees + double lambda = sunTrueLongitude - 0.00569 - 0.00478 * Math.sin(Math.toRadians(omega)); + return lambda; // in degrees } /** @@ -271,7 +253,6 @@ private static double getMeanObliquityOfEcliptic(double julianCenturies) { */ private static double getObliquityCorrection(double julianCenturies) { double obliquityOfEcliptic = getMeanObliquityOfEcliptic(julianCenturies); - double omega = 125.04 - 1934.136 * julianCenturies; return obliquityOfEcliptic + 0.00256 * Math.cos(Math.toRadians(omega)); // in degrees } @@ -288,9 +269,9 @@ private static double getObliquityCorrection(double julianCenturies) { private static double getSunDeclination(double julianCenturies) { double obliquityCorrection = getObliquityCorrection(julianCenturies); double lambda = getSunApparentLongitude(julianCenturies); - double sint = Math.sin(Math.toRadians(obliquityCorrection)) * Math.sin(Math.toRadians(lambda)); - return Math.toDegrees(Math.asin(sint)); // in degrees + double theta = Math.toDegrees(Math.asin(sint)); + return theta; // in degrees } /** @@ -307,62 +288,42 @@ private static double getEquationOfTime(double julianCenturies) { double geomMeanLongSun = getSunGeometricMeanLongitude(julianCenturies); double eccentricityEarthOrbit = getEarthOrbitEccentricity(julianCenturies); double geomMeanAnomalySun = getSunGeometricMeanAnomaly(julianCenturies); - double y = Math.tan(Math.toRadians(epsilon) / 2.0); y *= y; - double sin2l0 = Math.sin(2.0 * Math.toRadians(geomMeanLongSun)); double sinm = Math.sin(Math.toRadians(geomMeanAnomalySun)); double cos2l0 = Math.cos(2.0 * Math.toRadians(geomMeanLongSun)); double sin4l0 = Math.sin(4.0 * Math.toRadians(geomMeanLongSun)); double sin2m = Math.sin(2.0 * Math.toRadians(geomMeanAnomalySun)); - double equationOfTime = y * sin2l0 - 2.0 * eccentricityEarthOrbit * sinm + 4.0 * eccentricityEarthOrbit * y * sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * eccentricityEarthOrbit * eccentricityEarthOrbit * sin2m; return Math.toDegrees(equationOfTime) * 4.0; // in minutes of time } - + /** * Return the hour angle of the sun in - * radians at sunrise for the latitude. + * radians at for the latitude. * - * @param lat + * @param latitude * the latitude of observer in degrees - * @param solarDec + * @param solarDeclination * the declination angle of sun in degrees * @param zenith * the zenith + * @param solarEvent + * If the hour angle is for sunrise or sunset * @return hour angle of sunrise in radians */ - private static double getSunHourAngleAtSunrise(double lat, double solarDec, double zenith) { - double latRad = Math.toRadians(lat); - double sdRad = Math.toRadians(solarDec); - - return (Math.acos(Math.cos(Math.toRadians(zenith)) / (Math.cos(latRad) * Math.cos(sdRad)) - Math.tan(latRad) - * Math.tan(sdRad))); // in radians - } - - /** - * Returns the hour angle of the sun in radiansat sunset for the latitude. - * @todo use - {@link #getSunHourAngleAtSunrise(double, double, double)} implementation to avoid - * duplication of code. - * - * @param lat - * the latitude of observer in degrees - * @param solarDec - * the declination angle of sun in degrees - * @param zenith - * the zenith - * @return the hour angle of sunset in radians - */ - private static double getSunHourAngleAtSunset(double lat, double solarDec, double zenith) { - double latRad = Math.toRadians(lat); - double sdRad = Math.toRadians(solarDec); - + private static double getSunHourAngle(double latitude, double solarDeclination, double zenith, SolarEvent solarEvent) { + double latRad = Math.toRadians(latitude); + double sdRad = Math.toRadians(solarDeclination); double hourAngle = (Math.acos(Math.cos(Math.toRadians(zenith)) / (Math.cos(latRad) * Math.cos(sdRad)) - Math.tan(latRad) * Math.tan(sdRad))); - return -hourAngle; // in radians + + if (solarEvent == SolarEvent.SUNSET) { + hourAngle = -hourAngle; + } + return hourAngle; } /** @@ -370,29 +331,26 @@ private static double getSunHourAngleAtSunset(double lat, double solarDec, doubl * horizontal coordinate system at the given location at the given time. Can be negative if the sun is below the * horizon. Not corrected for altitude. * - * @param cal + * @param calendar * time of calculation - * @param lat + * @param latitude * latitude of location for calculation - * @param lon + * @param longitude * longitude of location for calculation * @return solar elevation in degrees - horizon is 0 degrees, civil twilight is -6 degrees */ - public static double getSolarElevation(Calendar cal, double lat, double lon) { - double julianDay = getJulianDay(cal); + public static double getSolarElevation(Calendar calendar, double latitude, double longitude) { + double julianDay = getJulianDay(calendar); double julianCenturies = getJulianCenturiesFromJulianDay(julianDay); - - double eot = getEquationOfTime(julianCenturies); - - double longitude = (cal.get(Calendar.HOUR_OF_DAY) + 12.0) - + (cal.get(Calendar.MINUTE) + eot + cal.get(Calendar.SECOND) / 60.0) / 60.0; - - longitude = -(longitude * 360.0 / 24.0) % 360.0; - double hourAngle_rad = Math.toRadians(lon - longitude); + Double eot = getEquationOfTime(julianCenturies); + double adjustedLongitude = (calendar.get(Calendar.HOUR_OF_DAY) + 12.0) + + (calendar.get(Calendar.MINUTE) + eot + calendar.get(Calendar.SECOND) / 60.0) / 60.0; + adjustedLongitude = -(adjustedLongitude * 360.0 / 24.0) % 360.0; + double hourAngle_rad = Math.toRadians(longitude - adjustedLongitude); double declination = getSunDeclination(julianCenturies); double dec_rad = Math.toRadians(declination); - double lat_rad = Math.toRadians(lat); + double lat_rad = Math.toRadians(latitude); return Math.toDegrees(Math.asin((Math.sin(lat_rad) * Math.sin(dec_rad)) + (Math.cos(lat_rad) * Math.cos(dec_rad) * Math.cos(hourAngle_rad)))); @@ -403,30 +361,26 @@ public static double getSolarElevation(Calendar cal, double lat, double lon) { * horizontal coordinate system at the given location at the given time. Not corrected for altitude. True south is 0 * degrees. * - * @param cal + * @param calendar * time of calculation - * @param lat + * @param latitude * latitude of location for calculation - * @param lon + * @param longitude * longitude of location for calculation * @return the solar azimuth */ - public static double getSolarAzimuth(Calendar cal, double lat, double lon) { - double julianDay = getJulianDay(cal); + public static double getSolarAzimuth(Calendar calendar, double latitude, double longitude) { + double julianDay = getJulianDay(calendar); double julianCenturies = getJulianCenturiesFromJulianDay(julianDay); - - double eot = getEquationOfTime(julianCenturies); - - double longitude = (cal.get(Calendar.HOUR_OF_DAY) + 12.0) - + (cal.get(Calendar.MINUTE) + eot + cal.get(Calendar.SECOND) / 60.0) / 60.0; - - longitude = -(longitude * 360.0 / 24.0) % 360.0; - double hourAngle_rad = Math.toRadians(lon - longitude); + Double eot = getEquationOfTime(julianCenturies); + double adjustedLongitude = (calendar.get(Calendar.HOUR_OF_DAY) + 12.0) + + (calendar.get(Calendar.MINUTE) + eot + calendar.get(Calendar.SECOND) / 60.0) / 60.0; + adjustedLongitude = -(adjustedLongitude * 360.0 / 24.0) % 360.0; + double hourAngle_rad = Math.toRadians(longitude - adjustedLongitude); double declination = getSunDeclination(julianCenturies); double dec_rad = Math.toRadians(declination); - double lat_rad = Math.toRadians(lat); - + double lat_rad = Math.toRadians(latitude); return Math.toDegrees(Math.atan(Math.sin(hourAngle_rad) / ((Math.cos(hourAngle_rad) * Math.sin(lat_rad)) - (Math.tan(dec_rad) * Math.cos(lat_rad)))))+180; @@ -451,30 +405,26 @@ private static double getSunriseUTC(double julianDay, double latitude, double lo // Find the time of solar noon at the location, and use that declination. This is better than start of the // Julian day - double noonmin = getSolarNoonUTC(julianCenturies, longitude); double tnoon = getJulianCenturiesFromJulianDay(julianDay + noonmin / 1440.0); // First pass to approximate sunrise (using solar noon) - - double eqTime = getEquationOfTime(tnoon); - double solarDec = getSunDeclination(tnoon); - double hourAngle = getSunHourAngleAtSunrise(latitude, solarDec, zenith); - + double equationOfTime = getEquationOfTime(tnoon); + double solarDeclination = getSunDeclination(tnoon); + double hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, SolarEvent.SUNRISE); double delta = longitude - Math.toDegrees(hourAngle); double timeDiff = 4 * delta; // in minutes of time - double timeUTC = 720 + timeDiff - eqTime; // in minutes + double timeUTC = 720 + timeDiff - equationOfTime; // in minutes // Second pass includes fractional Julian Day in gamma calc - double newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + timeUTC / 1440.0); - eqTime = getEquationOfTime(newt); - solarDec = getSunDeclination(newt); - hourAngle = getSunHourAngleAtSunrise(latitude, solarDec, zenith); + equationOfTime = getEquationOfTime(newt); + solarDeclination = getSunDeclination(newt); + hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, SolarEvent.SUNRISE); delta = longitude - Math.toDegrees(hourAngle); timeDiff = 4 * delta; - timeUTC = 720 + timeDiff - eqTime; // in minutes + timeUTC = 720 + timeDiff - equationOfTime; // in minutes return timeUTC; } @@ -498,23 +448,14 @@ private static double getSunriseUTC(double julianDay, double latitude, double lo public double getUTCNoon(Calendar calendar, GeoLocation geoLocation) { double julianDay = getJulianDay(calendar); double julianCenturies = getJulianCenturiesFromJulianDay(julianDay); - double noon = getSolarNoonUTC(julianCenturies, -geoLocation.getLongitude()); noon = noon / 60; - - // ensure that the time is >= 0 and < 24 - while (noon < 0.0) { - noon += 24.0; - } - while (noon >= 24.0) { - noon -= 24.0; - } - return noon; + return noon > 0 ? noon % 24 : noon % 24 + 24; // ensure that the time is >= 0 and < 24 } /** * Return the Universal Coordinated Time (UTC) - * of solar noon for the given day at the given location + * of of solar noon for the given day at the given location * on earth. * * @param julianCenturies @@ -532,14 +473,12 @@ private static double getSolarNoonUTC(double julianCenturies, double longitude) // First pass uses approximate solar noon to calculate equation of time double tnoon = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + longitude / 360.0); - double eqTime = getEquationOfTime(tnoon); - double solNoonUTC = 720 + (longitude * 4) - eqTime; // min - + double equationOfTime = getEquationOfTime(tnoon); + double solNoonUTC = 720 + (longitude * 4) - equationOfTime; // minutes double newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) - 0.5 + solNoonUTC / 1440.0); - - eqTime = getEquationOfTime(newt); - return 720 + (longitude * 4) - eqTime; // min + equationOfTime = getEquationOfTime(newt); + return 720 + (longitude * 4) - equationOfTime; // minutes } /** @@ -561,31 +500,26 @@ private static double getSunsetUTC(double julianDay, double latitude, double lon // Find the time of solar noon at the location, and use that declination. This is better than start of the // Julian day - double noonmin = getSolarNoonUTC(julianCenturies, longitude); double tnoon = getJulianCenturiesFromJulianDay(julianDay + noonmin / 1440.0); // First calculates sunrise and approx length of day - - double eqTime = getEquationOfTime(tnoon); - double solarDec = getSunDeclination(tnoon); - double hourAngle = getSunHourAngleAtSunset(latitude, solarDec, zenith); - + double equationOfTime = getEquationOfTime(tnoon); + double solarDeclination = getSunDeclination(tnoon); + double hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, SolarEvent.SUNSET); double delta = longitude - Math.toDegrees(hourAngle); double timeDiff = 4 * delta; - double timeUTC = 720 + timeDiff - eqTime; + double timeUTC = 720 + timeDiff - equationOfTime; // Second pass includes fractional Julian Day in gamma calc - double newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + timeUTC / 1440.0); - eqTime = getEquationOfTime(newt); - solarDec = getSunDeclination(newt); - hourAngle = getSunHourAngleAtSunset(latitude, solarDec, zenith); - + equationOfTime = getEquationOfTime(newt); + solarDeclination = getSunDeclination(newt); + hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, SolarEvent.SUNSET); delta = longitude - Math.toDegrees(hourAngle); timeDiff = 4 * delta; - timeUTC = 720 + timeDiff - eqTime; // in minutes + timeUTC = 720 + timeDiff - equationOfTime; // in minutes return timeUTC; } }