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 Jul 5, 2024
1 parent 44f2655 commit ba98e54
Show file tree
Hide file tree
Showing 4 changed files with 373 additions and 94 deletions.
4 changes: 3 additions & 1 deletion src/main/java/com/kosherjava/zmanim/AstronomicalCalendar.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,8 @@ public class AstronomicalCalendar implements Cloneable {
*/
private AstronomicalCalculator astronomicalCalculator;

private final TimeZone timeZoneUTC = TimeZone.getTimeZone("UTC");

/**
* The getSunrise method returns a <code>Date</code> representing the
* {@link AstronomicalCalculator#getElevationAdjustment(double) elevation adjusted} sunrise time. The zenith used
Expand Down Expand Up @@ -616,7 +618,7 @@ protected Date getDateFromTime(double time, SolarEvent solarEvent) {
double calculatedTime = time;

Calendar adjustedCalendar = getAdjustedCalendar();
Calendar cal = Calendar.getInstance(TimeZone.getTimeZone("UTC"));
Calendar cal = Calendar.getInstance(timeZoneUTC);
cal.clear();// clear all fields
cal.set(Calendar.YEAR, adjustedCalendar.get(Calendar.YEAR));
cal.set(Calendar.MONTH, adjustedCalendar.get(Calendar.MONTH));
Expand Down
245 changes: 173 additions & 72 deletions src/main/java/com/kosherjava/zmanim/util/NOAACalculator.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -311,11 +311,13 @@ private static double getEquationOfTime(double 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 geomMeanLongSunRad = Math.toRadians(geomMeanLongSun);
double geomMeanAnomalySunRad = Math.toRadians(geomMeanAnomalySun);
double sin2l0 = Math.sin(2.0 * geomMeanLongSunRad);
double sinm = Math.sin(geomMeanAnomalySunRad);
double cos2l0 = Math.cos(2.0 * geomMeanLongSunRad);
double sin4l0 = Math.sin(4.0 * geomMeanLongSunRad);
double sin2m = Math.sin(2.0 * geomMeanAnomalySunRad);

double equationOfTime = y * sin2l0 - 2.0 * eccentricityEarthOrbit * sinm + 4.0 * eccentricityEarthOrbit * y
* sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * eccentricityEarthOrbit * eccentricityEarthOrbit * sin2m;
Expand All @@ -337,34 +339,12 @@ private static double getEquationOfTime(double julianCenturies) {
private static double getSunHourAngleAtSunrise(double lat, double solarDec, double zenith) {
double latRad = Math.toRadians(lat);
double sdRad = Math.toRadians(solarDec);
double zRad = Math.toRadians(zenith);

return (Math.acos(Math.cos(Math.toRadians(zenith)) / (Math.cos(latRad) * Math.cos(sdRad)) - Math.tan(latRad)
return (Math.acos(Math.cos(zRad) / (Math.cos(latRad) * Math.cos(sdRad)) - Math.tan(latRad)
* Math.tan(sdRad))); // in radians
}

/**
* Returns the <a href="https://en.wikipedia.org/wiki/Hour_angle">hour angle</a> of the sun in <a href=
* "https://en.wikipedia.org/wiki/Radian">radians</a>at 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 <a href="https://en.wikipedia.org/wiki/Radian">radians</a>
*/
private static double getSunHourAngleAtSunset(double lat, double solarDec, double zenith) {
double latRad = Math.toRadians(lat);
double sdRad = Math.toRadians(solarDec);

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
}

/**
* Return the <a href="https://en.wikipedia.org/wiki/Celestial_coordinate_system">Solar Elevation</a> for the
* horizontal coordinate system at the given location at the given time. Can be negative if the sun is below the
Expand Down Expand Up @@ -447,35 +427,7 @@ public static double getSolarAzimuth(Calendar cal, double lat, double lon) {
* @return the time in minutes from zero UTC
*/
private static double getSunriseUTC(double julianDay, double latitude, double longitude, double zenith) {
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 eqTime = getEquationOfTime(tnoon);
double solarDec = getSunDeclination(tnoon);
double hourAngle = getSunHourAngleAtSunrise(latitude, solarDec, zenith);

double delta = longitude - Math.toDegrees(hourAngle);
double timeDiff = 4 * delta; // in minutes of time
double timeUTC = 720 + timeDiff - eqTime; // 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);
delta = longitude - Math.toDegrees(hourAngle);
timeDiff = 4 * delta;
timeUTC = 720 + timeDiff - eqTime; // in minutes
return timeUTC;
return getTimeUTC(julianDay, latitude, longitude, zenith, false);
}

/**
Expand Down Expand Up @@ -557,35 +509,184 @@ private static double getSolarNoonUTC(double julianCenturies, double longitude)
* @return the time in minutes from zero Universal Coordinated Time (UTC)
*/
private static double getSunsetUTC(double julianDay, double latitude, double longitude, double zenith) {
return getTimeUTC(julianDay, latitude, longitude, zenith, true);
}

/**
* Return the <a href="http://en.wikipedia.org/wiki/Universal_Coordinated_Time">Universal Coordinated Time</a> (UTC)
* of sunrise/sunset for the given day at the given location on earth
*
* @param julianDay
* the Julian day
* @param latitude
* the latitude of observer in degrees
* @param longitude
* the longitude of observer in degrees
* @param zenith
* the zenith
* @param isSunset
* True for sunset and false for sunrise.
* @return the time in minutes from zero UTC
*/
private static double getTimeUTC(double julianDay, double latitude, double longitude, double zenith, boolean isSunset) {
double julianCenturies = getJulianCenturiesFromJulianDay(julianDay);

// First pass to approximate sunrise (using solar noon)

double hourAngle = getHourAngleSunrise(julianDay, latitude, longitude, zenith);
if (Double.isNaN(hourAngle)) {
hourAngle = interpolateHourAngleSunrise(julianDay, latitude, longitude, zenith);
if (Double.isNaN(hourAngle)) return Double.NaN;
}
if (isSunset) 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);

// First calculates sunrise and approx length of day

double eqTime = getEquationOfTime(tnoon);
double solarDec = getSunDeclination(tnoon);
double hourAngle = getSunHourAngleAtSunset(latitude, solarDec, zenith);

double delta = longitude - Math.toDegrees(hourAngle);
double timeDiff = 4 * delta;
double timeUTC = 720 + timeDiff - eqTime;
double timeDiff = 4 * delta; // in minutes
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);
eqTime = getEquationOfTime(newt);
solarDec = getSunDeclination(newt);
hourAngle = getSunHourAngleAtSunset(latitude, solarDec, zenith);
hourAngle = getHourAngleSunset(julianDay, latitude, zenith, timeUTC);
if (Double.isNaN(hourAngle)) {
hourAngle = interpolateHourAngleSunset(julianDay, latitude, zenith, timeUTC);
if (Double.isNaN(hourAngle)) return Double.NaN;
}
if (isSunset) 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 interpolateHourAngleSunrise(double julianDay, double latitude, double longitude, double zenith) {
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 = getHourAngleSunrise(d, latitude, longitude, zenith);

if (!Double.isNaN(hourAngle)) {
x1 = d;
y1 = hourAngle;
break;
}
d--;
}

d = julianDay + 1;
final double dayLast = julianDay + 366;
while (d <= dayLast) {
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) {
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);
}
}
Loading

0 comments on commit ba98e54

Please sign in to comment.