Skip to content

Commit

Permalink
Interpolate missing dates
Browse files Browse the repository at this point in the history
  • Loading branch information
@pnemonic78
pnemonic78 committed Jun 25, 2022
1 parent 871c210 commit d7813f7
Show file tree
Hide file tree
Showing 2 changed files with 108 additions and 83 deletions.
141 changes: 90 additions & 51 deletions src/main/java/com/kosherjava/zmanim/util/NOAACalculator.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -480,76 +480,73 @@ 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);
double timeUTC = 720 + timeDiff - eqTime; // in minutes

// 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);
timeUTC = 720 + timeDiff - eqTime; // in minutes
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;
Expand All @@ -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) {
Expand All @@ -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);
}
}
50 changes: 18 additions & 32 deletions src/main/java/com/kosherjava/zmanim/util/SunTimesCalculator.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand Down

0 comments on commit d7813f7

Please sign in to comment.