diff --git a/src/main/java/com/kosherjava/zmanim/util/NOAACalculator.java b/src/main/java/com/kosherjava/zmanim/util/NOAACalculator.java index 3626e911..ba7cca2c 100644 --- a/src/main/java/com/kosherjava/zmanim/util/NOAACalculator.java +++ b/src/main/java/com/kosherjava/zmanim/util/NOAACalculator.java @@ -480,22 +480,20 @@ private static double getSunsetUTC(double julianDay, double latitude, double lon private static double getTimeUTC(double julianDay, double latitude, double longitude, double zenith, boolean isSunrise) { double julianCenturies = getJulianCenturiesFromJulianDay(julianDay); - // 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 solarDec = getSunDeclination(tnoon); - double hourAngle = getSunHourAngleAtSunrise(latitude, solarDec, zenith); + double hourAngle = getHourAngleSunrise(julianDay, latitude, longitude, zenith); if (Double.isNaN(hourAngle)) { - hourAngle = interpolateHourAngle(julianDay, latitude, longitude, zenith, Double.NaN); + hourAngle = interpolateHourAngleSunrise(julianDay, latitude, longitude, zenith); if (Double.isNaN(hourAngle)) return Double.NaN; } if (!isSunrise) hourAngle = -hourAngle; + // 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); double delta = longitude - Math.toDegrees(hourAngle); double timeDiff = 4 * delta; // in minutes of time double eqTime = getEquationOfTime(tnoon); @@ -503,16 +501,15 @@ private static double getTimeUTC(double julianDay, double latitude, double longi // Second pass includes fractional Julian Day in gamma calc - double newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + timeUTC - / 1440.0); - solarDec = getSunDeclination(newt); - hourAngle = getSunHourAngleAtSunrise(latitude, solarDec, zenith); + hourAngle = getHourAngleSunset(julianDay, latitude, zenith, timeUTC); if (Double.isNaN(hourAngle)) { - hourAngle = interpolateHourAngle(julianDay, latitude, longitude, zenith, timeUTC); + hourAngle = interpolateHourAngleSunset(julianDay, latitude, zenith, timeUTC); if (Double.isNaN(hourAngle)) return Double.NaN; } if (!isSunrise) hourAngle = -hourAngle; + double newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + timeUTC + / 1440.0); delta = longitude - Math.toDegrees(hourAngle); timeDiff = 4 * delta; eqTime = getEquationOfTime(newt); @@ -520,36 +517,36 @@ private static double getTimeUTC(double julianDay, double latitude, double longi return timeUTC; } + private static double getHourAngleSunrise(double julianDay, double latitude, double longitude, double zenith) { + double julianCenturies = getJulianCenturiesFromJulianDay(julianDay); + double noonmin = getSolarNoonUTC(julianCenturies, longitude); + double tnoon = getJulianCenturiesFromJulianDay(julianDay + noonmin / 1440.0); + double solarDec = getSunDeclination(tnoon); + return getSunHourAngleAtSunrise(latitude, solarDec, zenith); + } + + private static double getHourAngleSunset(double julianDay, double latitude, double zenith, double timeUTC) { + double julianCenturies = getJulianCenturiesFromJulianDay(julianDay); + double newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + timeUTC + / 1440.0); + double solarDec = getSunDeclination(newt); + return getSunHourAngleAtSunrise(latitude, solarDec, zenith); + } + /** * Use linear interpolation to calculate a missing angle. */ - private static double interpolateHourAngle(double julianDay, double latitude, double longitude, double zenith, double timeUTC) { - double dayOfYear = julianDay; + private static double interpolateHourAngleSunrise(double julianDay, double latitude, double longitude, double zenith) { + double hourAngle; double x1 = 0; double y1 = 0; double x2 = 0; double y2 = 0; - double julianCenturies; - double noonmin; - double tnoon; - double solarDec; - double hourAngle = 0; - double newt; - - double d = dayOfYear - 1; - while (d >= 1) { - julianCenturies = getJulianCenturiesFromJulianDay(d); - noonmin = getSolarNoonUTC(julianCenturies, longitude); - tnoon = getJulianCenturiesFromJulianDay(d + noonmin / 1440.0); - if (Double.isNaN(timeUTC)) { - solarDec = getSunDeclination(tnoon); - hourAngle = getSunHourAngleAtSunrise(latitude, solarDec, zenith); - } else { - newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + timeUTC - / 1440.0); - solarDec = getSunDeclination(newt); - hourAngle = getSunHourAngleAtSunrise(latitude, solarDec, zenith); - } + + double d = julianDay - 1; + final double dayFirst = Math.max(julianDay - 366, 1); + while (d >= dayFirst) { + hourAngle = getHourAngleSunrise(d, latitude, longitude, zenith); if (!Double.isNaN(hourAngle)) { x1 = d; @@ -559,21 +556,63 @@ private static double interpolateHourAngle(double julianDay, double latitude, do d--; } - d = dayOfYear + 1; - final double dayLast = dayOfYear + 366; + d = julianDay + 1; + final double dayLast = julianDay + 366; while (d <= dayLast) { - julianCenturies = getJulianCenturiesFromJulianDay(d); - noonmin = getSolarNoonUTC(julianCenturies, longitude); - tnoon = getJulianCenturiesFromJulianDay(d + noonmin / 1440.0); - if (Double.isNaN(timeUTC)) { - solarDec = getSunDeclination(tnoon); - hourAngle = getSunHourAngleAtSunrise(latitude, solarDec, zenith); - } else { - newt = getJulianCenturiesFromJulianDay(getJulianDayFromJulianCenturies(julianCenturies) + timeUTC - / 1440.0); - solarDec = getSunDeclination(newt); - hourAngle = getSunHourAngleAtSunrise(latitude, solarDec, zenith); + hourAngle = getHourAngleSunrise(d, latitude, longitude, zenith); + + if (!Double.isNaN(hourAngle)) { + if (x1 == 0) { + x1 = d; + y1 = hourAngle; + d++; + continue; + } + x2 = d; + y2 = hourAngle; + break; } + d++; + } + + if ((x1 == 0) || (x2 == 0)) { + return Double.NaN; + } + double dx = x2 - x1; + if (dx == 0) { + return Double.NaN; + } + double dy = y2 - y1; + return y1 + ((julianDay - x1) * dy / dx); + } + + /** + * Use linear interpolation to calculate a missing angle. + */ + private static double interpolateHourAngleSunset(double julianDay, double latitude, double zenith, double timeUTC) { + double hourAngle; + double x1 = 0; + double y1 = 0; + double x2 = 0; + double y2 = 0; + + double d = julianDay - 1; + final double dayFirst = Math.max(julianDay - 366, 1); + while (d >= dayFirst) { + hourAngle = getHourAngleSunset(d, latitude, zenith, timeUTC); + + if (!Double.isNaN(hourAngle)) { + x1 = d; + y1 = hourAngle; + break; + } + d--; + } + + d = julianDay + 1; + final double dayLast = julianDay + 366; + while (d <= dayLast) { + hourAngle = getHourAngleSunset(d, latitude, zenith, timeUTC); if (!Double.isNaN(hourAngle)) { if (x1 == 0) { @@ -597,6 +636,6 @@ private static double interpolateHourAngle(double julianDay, double latitude, do return Double.NaN; } double dy = y2 - y1; - return y1 + ((dayOfYear - x1) * dy / dx); + return y1 + ((julianDay - x1) * dy / dx); } } diff --git a/src/main/java/com/kosherjava/zmanim/util/SunTimesCalculator.java b/src/main/java/com/kosherjava/zmanim/util/SunTimesCalculator.java index 913086e8..d80a9e7a 100644 --- a/src/main/java/com/kosherjava/zmanim/util/SunTimesCalculator.java +++ b/src/main/java/com/kosherjava/zmanim/util/SunTimesCalculator.java @@ -226,26 +226,22 @@ private static double getLocalMeanTime(double localHour, double sunRightAscensio */ private static double getTimeUTC(Calendar calendar, GeoLocation geoLocation, double zenith, boolean isSunrise) { int dayOfYear = calendar.get(Calendar.DAY_OF_YEAR); - double sunMeanAnomaly = getMeanAnomaly(dayOfYear, geoLocation.getLongitude(), isSunrise); - double sunTrueLong = getSunTrueLongitude(sunMeanAnomaly); - double cosLocalHourAngle = getCosLocalHourAngle(sunTrueLong, geoLocation.getLatitude(), zenith); + double latitude = geoLocation.getLatitude(); + double longitude = geoLocation.getLongitude(); - double hourAngle; - if (isSunrise) { - hourAngle = 360.0 - acosDeg(cosLocalHourAngle); - } else { // sunset - hourAngle = acosDeg(cosLocalHourAngle); - } + double hourAngle = getHourAngle(dayOfYear, latitude, longitude, zenith, isSunrise); if (Double.isNaN(hourAngle)) { - hourAngle = interpolateHourAngle(calendar, geoLocation, zenith, isSunrise); + hourAngle = interpolateHourAngle(dayOfYear, geoLocation, zenith, isSunrise); if (Double.isNaN(hourAngle)) return Double.NaN; } double localHour = hourAngle / DEG_PER_HOUR; + double sunMeanAnomaly = getMeanAnomaly(dayOfYear, longitude, isSunrise); + double sunTrueLong = getSunTrueLongitude(sunMeanAnomaly); double sunRightAscensionHours = getSunRightAscensionHours(sunTrueLong); double localMeanTime = getLocalMeanTime(localHour, sunRightAscensionHours, - getApproxTimeDays(dayOfYear, getHoursFromMeridian(geoLocation.getLongitude()), isSunrise)); - double processedTime = localMeanTime - getHoursFromMeridian(geoLocation.getLongitude()); + getApproxTimeDays(dayOfYear, getHoursFromMeridian(longitude), isSunrise)); + double processedTime = localMeanTime - getHoursFromMeridian(longitude); while (processedTime < 0.0) { processedTime += 24.0; } @@ -255,11 +251,17 @@ private static double getTimeUTC(Calendar calendar, GeoLocation geoLocation, dou return processedTime; } + private static double getHourAngle(int dayOfYear, double latitude, double longitude, double zenith, boolean isSunrise) { + double sunMeanAnomaly = getMeanAnomaly(dayOfYear, longitude, isSunrise); + double sunTrueLong = getSunTrueLongitude(sunMeanAnomaly); + double cosLocalHourAngle = getCosLocalHourAngle(sunTrueLong, latitude, zenith); + return isSunrise ? (360.0 - acosDeg(cosLocalHourAngle)) : acosDeg(cosLocalHourAngle); + } + /** * Use linear interpolation to calculate a missing angle. */ - private static double interpolateHourAngle(Calendar calendar, GeoLocation geoLocation, double zenith, boolean isSunrise) { - int dayOfYear = calendar.get(Calendar.DAY_OF_YEAR); + private static double interpolateHourAngle(int dayOfYear, GeoLocation geoLocation, double zenith, boolean isSunrise) { double latitude = geoLocation.getLatitude(); double longitude = geoLocation.getLongitude(); double hourAngle; @@ -271,15 +273,7 @@ private static double interpolateHourAngle(Calendar calendar, GeoLocation geoLoc int d = dayOfYear - 1; while (d >= 1) { - double sunMeanAnomaly = getMeanAnomaly(d, longitude, isSunrise); - double sunTrueLong = getSunTrueLongitude(sunMeanAnomaly); - double cosLocalHourAngle = getCosLocalHourAngle(sunTrueLong, latitude, zenith); - - if (isSunrise) { - hourAngle = 360.0 - acosDeg(cosLocalHourAngle); - } else { // sunset - hourAngle = acosDeg(cosLocalHourAngle); - } + hourAngle = getHourAngle(d, latitude, longitude, zenith, isSunrise); if (!Double.isNaN(hourAngle)) { if (isLastDayOfYear && (x2 == 0)) { x2 = d; @@ -296,15 +290,7 @@ private static double interpolateHourAngle(Calendar calendar, GeoLocation geoLoc d = dayOfYear + 1; while (d <= 366) { - double sunMeanAnomaly = getMeanAnomaly(d, longitude, isSunrise); - double sunTrueLong = getSunTrueLongitude(sunMeanAnomaly); - double cosLocalHourAngle = getCosLocalHourAngle(sunTrueLong, latitude, zenith); - - if (isSunrise) { - hourAngle = 360.0 - acosDeg(cosLocalHourAngle); - } else { // sunset - hourAngle = acosDeg(cosLocalHourAngle); - } + hourAngle = getHourAngle(d, latitude, longitude, zenith, isSunrise); if (!Double.isNaN(hourAngle)) { if (x1 == 0) { x1 = d;