Skip to content

Commit a94d6f6

Browse files
pnemonic78pnemonic78
committed
Interpolate missing dates
1 parent e88a456 commit a94d6f6

File tree

6 files changed

+266
-41
lines changed

6 files changed

+266
-41
lines changed

‎src/main/java/com/kosherjava/zmanim/AstronomicalCalendar.java

+3-1
Original file line numberDiff line numberDiff line change
@@ -106,6 +106,8 @@ public class AstronomicalCalendar implements Cloneable {
106106
*/
107107
private AstronomicalCalculator astronomicalCalculator;
108108

109+
private final TimeZone timeZoneUTC = TimeZone.getTimeZone("UTC");
110+
109111
/**
110112
* The getSunrise method returns a <code>Date</code> representing the
111113
* {@link AstronomicalCalculator#getElevationAdjustment(double) elevation adjusted} sunrise time. The zenith used
@@ -625,7 +627,7 @@ protected Date getDateFromTime(double time, SolarEvent solarEvent) {
625627
double calculatedTime = time;
626628

627629
Calendar adjustedCalendar = getAdjustedCalendar();
628-
Calendar cal = Calendar.getInstance(TimeZone.getTimeZone("UTC"));
630+
Calendar cal = Calendar.getInstance(timeZoneUTC);
629631
cal.clear();// clear all fields
630632
cal.set(Calendar.YEAR, adjustedCalendar.get(Calendar.YEAR));
631633
cal.set(Calendar.MONTH, adjustedCalendar.get(Calendar.MONTH));

‎src/main/java/com/kosherjava/zmanim/util/GeoLocation.java

+10-7
Original file line numberDiff line numberDiff line change
@@ -88,6 +88,9 @@ public class GeoLocation implements Cloneable {
8888
/** constant for milliseconds in an hour (3,600,000) */
8989
private static final long HOUR_MILLIS = MINUTE_MILLIS * 60;
9090

91+
private static final double PI2 = Math.PI * 2;
92+
private static final double PI_4 = Math.PI / 4;
93+
9194
/**
9295
* Method to return the elevation in Meters <b>above</b> sea level.
9396
*
@@ -441,7 +444,7 @@ private double vincentyInverseFormula(GeoLocation location, int formula) {
441444
double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
442445

443446
double lambda = L;
444-
double lambdaP = 2 * Math.PI;
447+
double lambdaP = PI2;
445448
double iterLimit = 20;
446449
double sinLambda = 0;
447450
double cosLambda = 0;
@@ -510,10 +513,10 @@ private double vincentyInverseFormula(GeoLocation location, int formula) {
510513
*/
511514
public double getRhumbLineBearing(GeoLocation location) {
512515
double dLon = Math.toRadians(location.getLongitude() - getLongitude());
513-
double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4)
514-
/ Math.tan(Math.toRadians(getLatitude()) / 2 + Math.PI / 4));
516+
double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + PI_4)
517+
/ Math.tan(Math.toRadians(getLatitude()) / 2 + PI_4));
515518
if (Math.abs(dLon) > Math.PI)
516-
dLon = dLon > 0 ? -(2 * Math.PI - dLon) : (2 * Math.PI + dLon);
519+
dLon = dLon > 0 ? -(PI2 - dLon) : (PI2 + dLon);
517520
return Math.toDegrees(Math.atan2(dLon, dPhi));
518521
}
519522

@@ -529,16 +532,16 @@ public double getRhumbLineDistance(GeoLocation location) {
529532
double earthRadius = 6378137; // Earth's radius in meters (WGS-84)
530533
double dLat = Math.toRadians(location.getLatitude()) - Math.toRadians(getLatitude());
531534
double dLon = Math.abs(Math.toRadians(location.getLongitude()) - Math.toRadians(getLongitude()));
532-
double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4)
533-
/ Math.tan(Math.toRadians(getLatitude()) / 2 + Math.PI / 4));
535+
double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + PI_4)
536+
/ Math.tan(Math.toRadians(getLatitude()) / 2 + PI_4));
534537
double q = dLat / dPhi;
535538

536539
if (!(Math.abs(q) <= Double.MAX_VALUE)) {
537540
q = Math.cos(Math.toRadians(getLatitude()));
538541
}
539542
// if dLon over 180° take shorter rhumb across 180° meridian:
540543
if (dLon > Math.PI) {
541-
dLon = 2 * Math.PI - dLon;
544+
dLon = PI2 - dLon;
542545
}
543546
double d = Math.sqrt(dLat * dLat + q * q * dLon * dLon);
544547
return d * earthRadius;

‎src/main/java/com/kosherjava/zmanim/util/GeoLocationUtils.java

+10-7
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,9 @@ public class GeoLocationUtils {
4343
*/
4444
private static int FINAL_BEARING = 2;
4545

46+
private static final double PI2 = Math.PI * 2;
47+
private static final double PI_4 = Math.PI / 4;
48+
4649
/**
4750
* Calculate the <a href="http://en.wikipedia.org/wiki/Great_circle">geodesic</a> initial bearing between this Object and
4851
* a second Object passed to this method using <a href="http://en.wikipedia.org/wiki/Thaddeus_Vincenty">Thaddeus
@@ -127,7 +130,7 @@ private static double vincentyFormula(GeoLocation location, GeoLocation destinat
127130
double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
128131

129132
double lambda = L;
130-
double lambdaP = 2 * Math.PI;
133+
double lambdaP = PI2;
131134
double iterLimit = 20;
132135
double sinLambda = 0;
133136
double cosLambda = 0;
@@ -211,10 +214,10 @@ private static double vincentyFormula(GeoLocation location, GeoLocation destinat
211214
public static double getRhumbLineBearing(GeoLocation location, GeoLocation destination) {
212215
double dLon = Math.toRadians(destination.getLongitude() - location.getLongitude());
213216
double dPhi = Math.log(Math.tan(Math.toRadians(destination.getLatitude())
214-
/ 2 + Math.PI / 4)
215-
/ Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4));
217+
/ 2 + PI_4)
218+
/ Math.tan(Math.toRadians(location.getLatitude()) / 2 + PI_4));
216219
if (Math.abs(dLon) > Math.PI)
217-
dLon = dLon > 0 ? -(2 * Math.PI - dLon) : (2 * Math.PI + dLon);
220+
dLon = dLon > 0 ? -(PI2 - dLon) : (PI2 + dLon);
218221
return Math.toDegrees(Math.atan2(dLon, dPhi));
219222
}
220223

@@ -232,16 +235,16 @@ public static double getRhumbLineDistance(GeoLocation location, GeoLocation dest
232235
double earthRadius = 6378137; // Earth's radius in meters (WGS-84)
233236
double dLat = Math.toRadians(location.getLatitude()) - Math.toRadians(destination.getLatitude());
234237
double dLon = Math.abs(Math.toRadians(location.getLongitude()) - Math.toRadians(destination.getLongitude()));
235-
double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4)
236-
/ Math.tan(Math.toRadians(destination.getLatitude()) / 2 + Math.PI / 4));
238+
double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + PI_4)
239+
/ Math.tan(Math.toRadians(destination.getLatitude()) / 2 + PI_4));
237240
double q = dLat / dPhi;
238241

239242
if (!(Math.abs(q) <= Double.MAX_VALUE)) {
240243
q = Math.cos(Math.toRadians(destination.getLatitude()));
241244
}
242245
// if dLon over 180° take shorter rhumb across 180° meridian:
243246
if (dLon > Math.PI) {
244-
dLon = 2 * Math.PI - dLon;
247+
dLon = PI2 - dLon;
245248
}
246249
double d = Math.sqrt(dLat * dLat + q * q * dLon * dLon);
247250
return d * earthRadius;

‎src/main/java/com/kosherjava/zmanim/util/NOAACalculator.java

+157-7
Original file line numberDiff line numberDiff line change
@@ -284,11 +284,13 @@ private static double getEquationOfTime(double julianCenturies) {
284284
double geomMeanAnomalySun = getSunGeometricMeanAnomaly(julianCenturies);
285285
double y = Math.tan(Math.toRadians(epsilon) / 2.0);
286286
y *= y;
287-
double sin2l0 = Math.sin(2.0 * Math.toRadians(geomMeanLongSun));
288-
double sinm = Math.sin(Math.toRadians(geomMeanAnomalySun));
289-
double cos2l0 = Math.cos(2.0 * Math.toRadians(geomMeanLongSun));
290-
double sin4l0 = Math.sin(4.0 * Math.toRadians(geomMeanLongSun));
291-
double sin2m = Math.sin(2.0 * Math.toRadians(geomMeanAnomalySun));
287+
double geomMeanLongSunRad = Math.toRadians(geomMeanLongSun);
288+
double geomMeanAnomalySunRad = Math.toRadians(geomMeanAnomalySun);
289+
double sin2l0 = Math.sin(2.0 * geomMeanLongSunRad);
290+
double sinm = Math.sin(geomMeanAnomalySunRad);
291+
double cos2l0 = Math.cos(2.0 * geomMeanLongSunRad);
292+
double sin4l0 = Math.sin(4.0 * geomMeanLongSunRad);
293+
double sin2m = Math.sin(2.0 * geomMeanAnomalySunRad);
292294
double equationOfTime = y * sin2l0 - 2.0 * eccentricityEarthOrbit * sinm + 4.0 * eccentricityEarthOrbit * y
293295
* sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * eccentricityEarthOrbit * eccentricityEarthOrbit * sin2m;
294296
return Math.toDegrees(equationOfTime) * 4.0;
@@ -311,7 +313,8 @@ private static double getEquationOfTime(double julianCenturies) {
311313
private static double getSunHourAngle(double latitude, double solarDeclination, double zenith, SolarEvent solarEvent) {
312314
double latRad = Math.toRadians(latitude);
313315
double sdRad = Math.toRadians(solarDeclination);
314-
double hourAngle = (Math.acos(Math.cos(Math.toRadians(zenith)) / (Math.cos(latRad) * Math.cos(sdRad))
316+
double zRad = Math.toRadians(zenith);
317+
double hourAngle = (Math.acos(Math.cos(zRad) / (Math.cos(latRad) * Math.cos(sdRad))
315318
- Math.tan(latRad) * Math.tan(sdRad)));
316319

317320
if (solarEvent == SolarEvent.SUNSET) {
@@ -493,19 +496,166 @@ private static double getSunRiseSetUTC(Calendar calendar, double latitude, doubl
493496
double equationOfTime = getEquationOfTime(tnoon);
494497
double solarDeclination = getSunDeclination(tnoon);
495498
double hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
499+
if (Double.isNaN(hourAngle)) {
500+
hourAngle = interpolateHourAngle1(julianDay, latitude, longitude, zenith, solarEvent);
501+
if (Double.isNaN(hourAngle)) {
502+
return Double.NaN;
503+
}
504+
}
496505
double delta = longitude - Math.toDegrees(hourAngle);
497506
double timeDiff = 4 * delta;
498507
double timeUTC = 720 + timeDiff - equationOfTime;
499508

500509
// Second pass includes fractional Julian Day in gamma calc
501510
double newt = getJulianCenturiesFromJulianDay(julianDay + timeUTC / 1440.0);
502511
equationOfTime = getEquationOfTime(newt);
503-
512+
504513
solarDeclination = getSunDeclination(newt);
505514
hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
515+
if (Double.isNaN(hourAngle)) {
516+
hourAngle = interpolateHourAngle2(julianDay, latitude, longitude, zenith, solarEvent);
517+
if (Double.isNaN(hourAngle)) {
518+
return Double.NaN;
519+
}
520+
}
506521
delta = longitude - Math.toDegrees(hourAngle);
507522
timeDiff = 4 * delta;
508523
timeUTC = 720 + timeDiff - equationOfTime;
509524
return timeUTC;
510525
}
526+
527+
private static double getSunHourAngle1(double julianDay, double latitude, double longitude, double zenith, SolarEvent solarEvent) {
528+
double noonmin = getSolarNoonMidnightUTC(julianDay, longitude, SolarEvent.NOON);
529+
double tnoon = getJulianCenturiesFromJulianDay(julianDay + noonmin / 1440.0);
530+
double solarDeclination = getSunDeclination(tnoon);
531+
return getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
532+
}
533+
534+
private static double getSunHourAngle2(double julianDay, double latitude, double longitude, double zenith, SolarEvent solarEvent) {
535+
double noonmin = getSolarNoonMidnightUTC(julianDay, longitude, SolarEvent.NOON);
536+
double tnoon = getJulianCenturiesFromJulianDay(julianDay + noonmin / 1440.0);
537+
538+
// First calculates sunrise and approximate length of day
539+
double equationOfTime = getEquationOfTime(tnoon);
540+
double solarDeclination = getSunDeclination(tnoon);
541+
double hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
542+
if (Double.isNaN(hourAngle)) {
543+
return Double.NaN;
544+
}
545+
double delta = longitude - Math.toDegrees(hourAngle);
546+
double timeDiff = 4 * delta;
547+
double timeUTC = 720 + timeDiff - equationOfTime;
548+
549+
// Second pass includes fractional Julian Day in gamma calc
550+
double newt = getJulianCenturiesFromJulianDay(julianDay + timeUTC / 1440.0);
551+
552+
solarDeclination = getSunDeclination(newt);
553+
return getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
554+
}
555+
556+
/**
557+
* Use linear interpolation to calculate a missing angle.
558+
*/
559+
private static double interpolateHourAngle1(double julianDay, double latitude, double longitude, double zenith, SolarEvent solarEvent) {
560+
double hourAngle;
561+
double x1 = 0;
562+
double y1 = 0;
563+
double x2 = 0;
564+
double y2 = 0;
565+
566+
double d = julianDay - 1;
567+
final double dayFirst = Math.max(julianDay - 366, 1);
568+
while (d >= dayFirst) {
569+
hourAngle = getSunHourAngle1(d, latitude, longitude, zenith, solarEvent);
570+
571+
if (!Double.isNaN(hourAngle)) {
572+
x1 = d;
573+
y1 = hourAngle;
574+
break;
575+
}
576+
d--;
577+
}
578+
579+
d = julianDay + 1;
580+
final double dayLast = julianDay + 366;
581+
while (d <= dayLast) {
582+
hourAngle = getSunHourAngle1(d, latitude, longitude, zenith, solarEvent);
583+
584+
if (!Double.isNaN(hourAngle)) {
585+
if (x1 == 0) {
586+
x1 = d;
587+
y1 = hourAngle;
588+
d++;
589+
continue;
590+
}
591+
x2 = d;
592+
y2 = hourAngle;
593+
break;
594+
}
595+
d++;
596+
}
597+
598+
if ((x1 == 0) || (x2 == 0)) {
599+
return Double.NaN;
600+
}
601+
double dx = x2 - x1;
602+
if (dx == 0) {
603+
return Double.NaN;
604+
}
605+
double dy = y2 - y1;
606+
return y1 + ((julianDay - x1) * dy / dx);
607+
}
608+
609+
/**
610+
* Use linear interpolation to calculate a missing angle.
611+
*/
612+
private static double interpolateHourAngle2(double julianDay, double latitude, double longitude, double zenith, SolarEvent solarEvent) {
613+
double hourAngle;
614+
double x1 = 0;
615+
double y1 = 0;
616+
double x2 = 0;
617+
double y2 = 0;
618+
619+
double d = julianDay - 1;
620+
final double dayFirst = Math.max(julianDay - 366, 1);
621+
while (d >= dayFirst) {
622+
hourAngle = getSunHourAngle2(d, latitude, longitude, zenith, solarEvent);
623+
624+
if (!Double.isNaN(hourAngle)) {
625+
x1 = d;
626+
y1 = hourAngle;
627+
break;
628+
}
629+
d--;
630+
}
631+
632+
d = julianDay + 1;
633+
final double dayLast = julianDay + 366;
634+
while (d <= dayLast) {
635+
hourAngle = getSunHourAngle2(d, latitude, longitude, zenith, solarEvent);
636+
637+
if (!Double.isNaN(hourAngle)) {
638+
if (x1 == 0) {
639+
x1 = d;
640+
y1 = hourAngle;
641+
d++;
642+
continue;
643+
}
644+
x2 = d;
645+
y2 = hourAngle;
646+
break;
647+
}
648+
d++;
649+
}
650+
651+
if ((x1 == 0) || (x2 == 0)) {
652+
return Double.NaN;
653+
}
654+
double dx = x2 - x1;
655+
if (dx == 0) {
656+
return Double.NaN;
657+
}
658+
double dy = y2 - y1;
659+
return y1 + ((julianDay - x1) * dy / dx);
660+
}
511661
}

0 commit comments

Comments
 (0)