Skip to content

Commit fd322f5

Browse files
pnemonic78pnemonic78
committed
Interpolate missing dates
1 parent 2f69b8c commit fd322f5

File tree

6 files changed

+273
-46
lines changed

6 files changed

+273
-46
lines changed

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

Lines changed: 3 additions & 1 deletion
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
@@ -616,7 +618,7 @@ protected Date getDateFromTime(double time, SolarEvent solarEvent) {
616618
double calculatedTime = time;
617619

618620
Calendar adjustedCalendar = getAdjustedCalendar();
619-
Calendar cal = Calendar.getInstance(TimeZone.getTimeZone("UTC"));
621+
Calendar cal = Calendar.getInstance(timeZoneUTC);
620622
cal.clear();// clear all fields
621623
cal.set(Calendar.YEAR, adjustedCalendar.get(Calendar.YEAR));
622624
cal.set(Calendar.MONTH, adjustedCalendar.get(Calendar.MONTH));

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

Lines changed: 10 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -83,6 +83,9 @@ public class GeoLocation implements Cloneable {
8383
/** constant for milliseconds in an hour (3,600,000) */
8484
private static final long HOUR_MILLIS = MINUTE_MILLIS * 60;
8585

86+
private static final double PI2 = Math.PI * 2;
87+
private static final double PI_4 = Math.PI / 4;
88+
8689
/**
8790
* Method to get the elevation in Meters.
8891
*
@@ -426,7 +429,7 @@ private double vincentyFormula(GeoLocation location, int formula) {
426429
double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
427430

428431
double lambda = L;
429-
double lambdaP = 2 * Math.PI;
432+
double lambdaP = PI2;
430433
double iterLimit = 20;
431434
double sinLambda = 0;
432435
double cosLambda = 0;
@@ -495,10 +498,10 @@ private double vincentyFormula(GeoLocation location, int formula) {
495498
*/
496499
public double getRhumbLineBearing(GeoLocation location) {
497500
double dLon = Math.toRadians(location.getLongitude() - getLongitude());
498-
double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4)
499-
/ Math.tan(Math.toRadians(getLatitude()) / 2 + Math.PI / 4));
501+
double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + PI_4)
502+
/ Math.tan(Math.toRadians(getLatitude()) / 2 + PI_4));
500503
if (Math.abs(dLon) > Math.PI)
501-
dLon = dLon > 0 ? -(2 * Math.PI - dLon) : (2 * Math.PI + dLon);
504+
dLon = dLon > 0 ? -(PI2 - dLon) : (PI2 + dLon);
502505
return Math.toDegrees(Math.atan2(dLon, dPhi));
503506
}
504507

@@ -514,16 +517,16 @@ public double getRhumbLineDistance(GeoLocation location) {
514517
double earthRadius = 6378137; // Earth's radius in meters (WGS-84)
515518
double dLat = Math.toRadians(location.getLatitude()) - Math.toRadians(getLatitude());
516519
double dLon = Math.abs(Math.toRadians(location.getLongitude()) - Math.toRadians(getLongitude()));
517-
double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + Math.PI / 4)
518-
/ Math.tan(Math.toRadians(getLatitude()) / 2 + Math.PI / 4));
520+
double dPhi = Math.log(Math.tan(Math.toRadians(location.getLatitude()) / 2 + PI_4)
521+
/ Math.tan(Math.toRadians(getLatitude()) / 2 + PI_4));
519522
double q = dLat / dPhi;
520523

521524
if (!(Math.abs(q) <= Double.MAX_VALUE)) {
522525
q = Math.cos(Math.toRadians(getLatitude()));
523526
}
524527
// if dLon over 180° take shorter rhumb across 180° meridian:
525528
if (dLon > Math.PI) {
526-
dLon = 2 * Math.PI - dLon;
529+
dLon = PI2 - dLon;
527530
}
528531
double d = Math.sqrt(dLat * dLat + q * q * dLon * dLon);
529532
return d * earthRadius;

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

Lines changed: 10 additions & 7 deletions
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‎

Lines changed: 163 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -276,11 +276,13 @@ private static double getEquationOfTime(double julianCenturies) {
276276
double geomMeanAnomalySun = getSunGeometricMeanAnomaly(julianCenturies);
277277
double y = Math.tan(Math.toRadians(epsilon) / 2.0);
278278
y *= y;
279-
double sin2l0 = Math.sin(2.0 * Math.toRadians(geomMeanLongSun));
280-
double sinm = Math.sin(Math.toRadians(geomMeanAnomalySun));
281-
double cos2l0 = Math.cos(2.0 * Math.toRadians(geomMeanLongSun));
282-
double sin4l0 = Math.sin(4.0 * Math.toRadians(geomMeanLongSun));
283-
double sin2m = Math.sin(2.0 * Math.toRadians(geomMeanAnomalySun));
279+
double geomMeanLongSunRad = Math.toRadians(geomMeanLongSun);
280+
double geomMeanAnomalySunRad = Math.toRadians(geomMeanAnomalySun);
281+
double sin2l0 = Math.sin(2.0 * geomMeanLongSunRad);
282+
double sinm = Math.sin(geomMeanAnomalySunRad);
283+
double cos2l0 = Math.cos(2.0 * geomMeanLongSunRad);
284+
double sin4l0 = Math.sin(4.0 * geomMeanLongSunRad);
285+
double sin2m = Math.sin(2.0 * geomMeanAnomalySunRad);
284286
double equationOfTime = y * sin2l0 - 2.0 * eccentricityEarthOrbit * sinm + 4.0 * eccentricityEarthOrbit * y
285287
* sinm * cos2l0 - 0.5 * y * y * sin4l0 - 1.25 * eccentricityEarthOrbit * eccentricityEarthOrbit * sin2m;
286288
return Math.toDegrees(equationOfTime) * 4.0;
@@ -303,7 +305,8 @@ private static double getEquationOfTime(double julianCenturies) {
303305
private static double getSunHourAngle(double latitude, double solarDeclination, double zenith, SolarEvent solarEvent) {
304306
double latRad = Math.toRadians(latitude);
305307
double sdRad = Math.toRadians(solarDeclination);
306-
double hourAngle = (Math.acos(Math.cos(Math.toRadians(zenith)) / (Math.cos(latRad) * Math.cos(sdRad))
308+
double zRad = Math.toRadians(zenith);
309+
double hourAngle = (Math.acos(Math.cos(zRad) / (Math.cos(latRad) * Math.cos(sdRad))
307310
- Math.tan(latRad) * Math.tan(sdRad)));
308311

309312
if (solarEvent == SolarEvent.SUNSET) {
@@ -453,21 +456,170 @@ private static double getSunRiseSetUTC(Calendar calendar, double latitude, doubl
453456

454457
// First calculates sunrise and approximate length of day
455458
double equationOfTime = getEquationOfTime(tnoon);
456-
double solarDeclination = getSunDeclination(tnoon);
457-
double hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
459+
double hourAngle = getSunHourAngle1(julianDay, latitude, longitude, zenith, solarEvent);
460+
if (Double.isNaN(hourAngle)) {
461+
hourAngle = interpolateHourAngle1(julianDay, latitude, longitude, zenith, solarEvent);
462+
if (Double.isNaN(hourAngle)) {
463+
return Double.NaN;
464+
}
465+
}
458466
double delta = longitude - Math.toDegrees(hourAngle);
459467
double timeDiff = 4 * delta;
460468
double timeUTC = 720 + timeDiff - equationOfTime;
461469

462470
// Second pass includes fractional Julian Day in gamma calc
463471
double newt = getJulianCenturiesFromJulianDay(julianDay + timeUTC / 1440.0);
464472
equationOfTime = getEquationOfTime(newt);
465-
466-
solarDeclination = getSunDeclination(newt);
467-
hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
473+
474+
hourAngle = getSunHourAngle2(julianDay, latitude, longitude, zenith, solarEvent);
475+
if (Double.isNaN(hourAngle)) {
476+
hourAngle = interpolateHourAngle2(julianDay, latitude, longitude, zenith, solarEvent);
477+
if (Double.isNaN(hourAngle)) {
478+
return Double.NaN;
479+
}
480+
}
468481
delta = longitude - Math.toDegrees(hourAngle);
469482
timeDiff = 4 * delta;
470483
timeUTC = 720 + timeDiff - equationOfTime;
471484
return timeUTC;
472485
}
486+
487+
private static double getSunHourAngle1(double julianDay, double latitude, double longitude, double zenith, SolarEvent solarEvent) {
488+
// Find the time of solar noon at the location, and use that declination.
489+
// This is better than start of the Julian day
490+
double noonmin = getSolarNoonUTC(julianDay, longitude);
491+
double tnoon = getJulianCenturiesFromJulianDay(julianDay + noonmin / 1440.0);
492+
493+
// First calculates sunrise and approximate length of day
494+
double solarDeclination = getSunDeclination(tnoon);
495+
496+
return getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
497+
}
498+
499+
private static double getSunHourAngle2(double julianDay, double latitude, double longitude, double zenith, SolarEvent solarEvent) {
500+
// Find the time of solar noon at the location, and use that declination.
501+
// This is better than start of the Julian day
502+
double noonmin = getSolarNoonUTC(julianDay, longitude);
503+
double tnoon = getJulianCenturiesFromJulianDay(julianDay + noonmin / 1440.0);
504+
505+
// First calculates sunrise and approximate length of day
506+
double equationOfTime = getEquationOfTime(tnoon);
507+
double solarDeclination = getSunDeclination(tnoon);
508+
double hourAngle = getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
509+
double delta = longitude - Math.toDegrees(hourAngle);
510+
double timeDiff = 4 * delta;
511+
double timeUTC = 720 + timeDiff - equationOfTime;
512+
513+
// Second pass includes fractional Julian Day in gamma calc
514+
double newt = getJulianCenturiesFromJulianDay(julianDay + timeUTC / 1440.0);
515+
516+
solarDeclination = getSunDeclination(newt);
517+
return getSunHourAngle(latitude, solarDeclination, zenith, solarEvent);
518+
}
519+
520+
/**
521+
* Use linear interpolation to calculate a missing angle.
522+
*/
523+
private static double interpolateHourAngle1(double julianDay, double latitude, double longitude, double zenith, SolarEvent solarEvent) {
524+
double hourAngle;
525+
double x1 = 0;
526+
double y1 = 0;
527+
double x2 = 0;
528+
double y2 = 0;
529+
530+
double d = julianDay - 1;
531+
final double dayFirst = Math.max(julianDay - 366, 1);
532+
while (d >= dayFirst) {
533+
hourAngle = getSunHourAngle1(d, latitude, longitude, zenith, solarEvent);
534+
535+
if (!Double.isNaN(hourAngle)) {
536+
x1 = d;
537+
y1 = hourAngle;
538+
break;
539+
}
540+
d--;
541+
}
542+
543+
d = julianDay + 1;
544+
final double dayLast = julianDay + 366;
545+
while (d <= dayLast) {
546+
hourAngle = getSunHourAngle1(d, latitude, longitude, zenith, solarEvent);
547+
548+
if (!Double.isNaN(hourAngle)) {
549+
if (x1 == 0) {
550+
x1 = d;
551+
y1 = hourAngle;
552+
d++;
553+
continue;
554+
}
555+
x2 = d;
556+
y2 = hourAngle;
557+
break;
558+
}
559+
d++;
560+
}
561+
562+
if ((x1 == 0) || (x2 == 0)) {
563+
return Double.NaN;
564+
}
565+
double dx = x2 - x1;
566+
if (dx == 0) {
567+
return Double.NaN;
568+
}
569+
double dy = y2 - y1;
570+
return y1 + ((julianDay - x1) * dy / dx);
571+
}
572+
573+
/**
574+
* Use linear interpolation to calculate a missing angle.
575+
*/
576+
private static double interpolateHourAngle2(double julianDay, double latitude, double longitude, double zenith, SolarEvent solarEvent) {
577+
double hourAngle;
578+
double x1 = 0;
579+
double y1 = 0;
580+
double x2 = 0;
581+
double y2 = 0;
582+
583+
double d = julianDay - 1;
584+
final double dayFirst = Math.max(julianDay - 366, 1);
585+
while (d >= dayFirst) {
586+
hourAngle = getSunHourAngle2(d, latitude, longitude, zenith, solarEvent);
587+
588+
if (!Double.isNaN(hourAngle)) {
589+
x1 = d;
590+
y1 = hourAngle;
591+
break;
592+
}
593+
d--;
594+
}
595+
596+
d = julianDay + 1;
597+
final double dayLast = julianDay + 366;
598+
while (d <= dayLast) {
599+
hourAngle = getSunHourAngle2(d, latitude, longitude, zenith, solarEvent);
600+
601+
if (!Double.isNaN(hourAngle)) {
602+
if (x1 == 0) {
603+
x1 = d;
604+
y1 = hourAngle;
605+
d++;
606+
continue;
607+
}
608+
x2 = d;
609+
y2 = hourAngle;
610+
break;
611+
}
612+
d++;
613+
}
614+
615+
if ((x1 == 0) || (x2 == 0)) {
616+
return Double.NaN;
617+
}
618+
double dx = x2 - x1;
619+
if (dx == 0) {
620+
return Double.NaN;
621+
}
622+
double dy = y2 - y1;
623+
return y1 + ((julianDay - x1) * dy / dx);
624+
}
473625
}

0 commit comments

Comments
 (0)