Skip to content

Commit

Permalink
astronomical sidereal longitude
Browse files Browse the repository at this point in the history 
see issue #851
  • Loading branch information
@MenoData
MenoData committed Jun 30, 2020
1 parent 7462703 commit 451120f
Show file tree
Hide file tree
Showing 4 changed files with 105 additions and 18 deletions.
16 changes: 8 additions & 8 deletions base/src/main/java/net/time4j/calendar/astro/StdSolarCalculator.java
View file
Expand Up @@ -496,9 +496,9 @@ private double equationOfCenter(double jct) {
* <ul>
* <li>right-ascension</li>
* <li>declination</li>
* <li>mean-anomaly</li>
* <li>nutation</li>
* <li>obliquity</li>
* <li>mean-anomaly</li>
* <li>solar-longitude</li>
* </ul>
*/
Expand All @@ -514,9 +514,9 @@ private double equationOfCenter(double jct) {
* <ul>
* <li>right-ascension</li>
* <li>declination</li>
* <li>mean-anomaly</li>
* <li>nutation</li>
* <li>obliquity</li>
* <li>mean-anomaly</li>
* <li>solar-longitude</li>
* </ul>
*/
Expand Down Expand Up @@ -593,12 +593,12 @@ public double getFeature(
double y = Math.cos(Math.toRadians(obliquity(jct))) * Math.sin(lRad);
double ra = Math.toDegrees(Math.atan2(y, Math.cos(lRad)));
return AstroUtils.toRange_0_360(ra);
case "mean-anomaly":
return meanAnomaly(jct);
case "nutation":
return nutation(jct);
case "obliquity":
return obliquity(jct);
case "mean-anomaly":
return meanAnomaly(jct);
case "solar-longitude":
return apparentSolarLongitude(jct, nutation(jct));
case "solar-latitude":
Expand Down Expand Up @@ -700,9 +700,9 @@ private double excentricity(double jct) {
* <ul>
* <li>right-ascension</li>
* <li>declination</li>
* <li>mean-anomaly</li>
* <li>nutation</li>
* <li>obliquity</li>
* <li>mean-anomaly</li>
* <li>solar-longitude</li>
* </ul>
*/
Expand All @@ -724,9 +724,9 @@ private double excentricity(double jct) {
* <ul>
* <li>right-ascension</li>
* <li>declination</li>
* <li>mean-anomaly</li>
* <li>nutation</li>
* <li>obliquity</li>
* <li>mean-anomaly</li>
* <li>solar-longitude</li>
* </ul>
*/
Expand Down Expand Up @@ -783,6 +783,8 @@ public double getFeature(
double ra = Math.toDegrees(Math.atan2(y, Math.cos(lRad)));
return AstroUtils.toRange_0_360(ra);
}
case "mean-anomaly":
return meanAnomaly(jct);
case "nutation": {
double[] result = new double[2];
nutations(jct, result);
Expand All @@ -793,8 +795,6 @@ public double getFeature(
nutations(jct, result);
return meanObliquity(jct) + result[1];
}
case "mean-anomaly":
return meanAnomaly(jct);
case "solar-longitude": {
double[] result = new double[2];
nutations(jct, result);
Expand Down
40 changes: 32 additions & 8 deletions base/src/main/java/net/time4j/calendar/hindu/HinduVariant.java
View file
Expand Up @@ -25,10 +25,11 @@
import net.time4j.PlainDate;
import net.time4j.calendar.IndianMonth;
import net.time4j.calendar.astro.GeoLocation;
import net.time4j.calendar.astro.JulianDay;
import net.time4j.calendar.astro.StdSolarCalculator;
import net.time4j.engine.CalendarDays;
import net.time4j.engine.EpochDays;
import net.time4j.engine.VariantSource;
import net.time4j.scale.TimeScale;

import java.io.IOException;
import java.io.InvalidObjectException;
Expand Down Expand Up @@ -773,12 +774,14 @@ private void readObject(ObjectInputStream in) throws IOException {
//~ Innere Klassen ----------------------------------------------------

@SuppressWarnings("ConstantConditions")
private static class ModernHinduCS
static class ModernHinduCS
extends HinduCS {

//~ Statische Felder/Initialisierungen ----------------------------

private static final double SIDEREAL_YEAR = 365d + (279457.0 / 1080000.0);
static final double SIDEREAL_YEAR = 365d + (279457.0 / 1080000.0);
static final double SIDEREAL_START = 336.13606783930777;

private static final double SIDEREAL_MONTH = 27d + (4644439.0 / 14438334.0);
private static final double SYNODIC_MONTH = 29d + (7087771.0 / 13358334.0);
private static final double EPSILON = Math.pow(2, -1000); // max recursion depth: 1000
Expand Down Expand Up @@ -995,10 +998,31 @@ private static double hTruePosition(
return modulo(lambda - equation, 360d);
}

private static double hSolarLongitude(double t) {
private static double hSiderealLongitude(double t) {
return modulo(hSolarLongitude(t) - hPrecession(t) + SIDEREAL_START, 360);
}

static double hSolarLongitude(double t) {
return hTruePosition(t, SIDEREAL_YEAR, 14d / 360d, ANOMALISTIC_YEAR, 1d / 42d);
}

// verified with Meeus example 21.c
static double hPrecession(double t) {
long unix = (long) (t + 1721424L - 2440587L) * 86400;
double jct = JulianDay.ofEphemerisTime(Moment.of(unix, TimeScale.POSIX)).getCenturyJ2000();

// using Meeus (21.6)
double eta = modulo(((47.0029 / 3600) + ((-0.03302 / 3600) + (0.00006 / 3600) * jct) * jct) * jct, 360);
double P = modulo(174.876384 + ((-869.8089 / 3600) + (0.03536 / 3600) * jct) * jct, 360);
double p = modulo(((5029.0966 / 3600) + ((1.11113 / 3600) + (-0.000006 / 3600) * jct) * jct) * jct, 360);

// use solar latitude = 0 and solar longitude (at mesha samkranti) = 0 in Meeus (21.7)
double A = Math.cos(Math.toRadians(eta)) * Math.sin(Math.toRadians(P));
double B = Math.cos(Math.toRadians(P));
double arg = Math.toDegrees(Math.atan2(A, B));
return modulo(p + P - arg, 360);
}

private static int hZodiac(double t) {
return (int) (Math.floor(hSolarLongitude(t) / 30d) + 1);
}
Expand All @@ -1017,10 +1041,10 @@ private static int hLunarDayFromMoment(double t) {

private static double hNewMoonBefore(double t) {
double tau = t - ((hLunarPhase(t) * SYNODIC_MONTH) / 360d);
return binarySearch(tau - 1, Math.min(t, tau + 1));
return binarySearchLunarPhase(tau - 1, Math.min(t, tau + 1));
}

private static double binarySearch(
private static double binarySearchLunarPhase(
double low,
double high
) {
Expand All @@ -1031,9 +1055,9 @@ private static double binarySearch(
}

if (hLunarPhase((low + high) / 2) < 180) {
return binarySearch(low, x);
return binarySearchLunarPhase(low, x);
} else {
return binarySearch(x, high);
return binarySearchLunarPhase(x, high);
}
}

Expand Down
65 changes: 64 additions & 1 deletion base/src/test/java/net/time4j/calendar/hindu/HinduVariantTest.java
View file
@@ -1,20 +1,29 @@
package net.time4j.calendar.hindu;

import net.time4j.Moment;
import net.time4j.PlainDate;
import net.time4j.PlainTimestamp;
import net.time4j.Weekday;
import net.time4j.calendar.IndianMonth;
import net.time4j.calendar.astro.GeoLocation;
import net.time4j.calendar.astro.JulianDay;
import net.time4j.calendar.astro.SolarTime;
import net.time4j.calendar.astro.StdSolarCalculator;
import net.time4j.engine.CalendarSystem;
import net.time4j.engine.EpochDays;
import net.time4j.scale.TimeScale;
import net.time4j.tz.ZonalOffset;
import org.junit.Test;
import org.junit.runner.RunWith;
import org.junit.runners.JUnit4;

import java.math.BigDecimal;
import java.util.Locale;
import java.util.concurrent.TimeUnit;

import static org.hamcrest.CoreMatchers.is;
import static org.junit.Assert.assertThat;
import static net.time4j.calendar.hindu.HinduVariant.ModernHinduCS.*;


@RunWith(JUnit4.class)
Expand Down Expand Up @@ -331,7 +340,15 @@ public void chennai() {

System.out.println(cal.transform(PlainDate.axis()));

int[] lengthOfMonths = {30, 32, 31, 32, 31, 30, 30, 30, 29, 30, 29, 31};
int[][] lengthOfMonths = {
{30, 32, 31, 32, 31, 30, 30, 30, 29, 30, 29, 31}, // 2020
{31, 31, 32, 31, 31, 31, 30, 29, 29, 30, 30, 30},
{31, 31, 32, 31, 31, 31, 30, 29, 30, 29, 30, 30},
{31, 31, 32, 31, 31, 31, 30, 29, 30, 29, 30, 31},
{30, 32, 31, 32, 31, 30, 30, 30, 29, 30, 30, 30},
// {31, 31, 31, 32, 31, 30, 30, 30, 29, }, // 2025
};

// 31, 32, 31, 32, 31, 30, 30, 30, 29, 30, 29, 31,
// tamil: 31, 31, 32, 31, 31, 31, 30, 29, 29, 30, 30, 30
// tamil-ujjain: 31, 31, 32, 31, 31, 31, 30, 29, 29, 30, 30, 30
Expand Down Expand Up @@ -369,4 +386,50 @@ public void chennai() {
}
}

@Test
public void meshaSamkranti285() { // CC, 20.41
long offset = ZonalOffset.atLongitude(new BigDecimal(HinduVariant.UJJAIN.getLongitude())).getIntegralAmount();
double fixed = ((meshaSamkranti(285) * 86400d) - offset) / 86400d;
Moment m = Moment.of(86400 * ((long) fixed + 1721424L - 2440587L), TimeScale.POSIX);
assertThat(
m.get(PlainDate.YEAR.atUTC()),
is(285));
assertThat(
Math.abs(SIDEREAL_START - hPrecession(fixed)) < 0.001,
is(true));
}

@Test
public void meshaSamkranti2000() { // CC-example, page 364
long fixed = (long) (meshaSamkranti(2000) * 86400d);
long unix = fixed + (1721424L - 2440587L) * 86400;
Moment s = Moment.of(unix, TimeScale.POSIX).with(Moment.PRECISION, TimeUnit.MINUTES);
assertThat(s, is(PlainTimestamp.of(2000, 4, 13, 17, 55).atUTC()));
}

private static double meshaSamkranti(int gyear) {
long t = PlainDate.of(gyear, 1, 1).get(EpochDays.RATA_DIE);
double tau = t + ((SIDEREAL_YEAR / 360d) * HinduCS.modulo(-hSolarLongitude(t), 360));
double a = Math.max(t, tau - 5);
double b = tau + 5;
return binarySearchSolarLongitude(a, b);
}

private static double binarySearchSolarLongitude(
double low,
double high
) {
double x = (low + high) / 2;

if (high - low < 0.00001) {
return x;
}

if (HinduCS.modulo(hSolarLongitude((low + high) / 2), 360) < 180) {
return binarySearchSolarLongitude(low, x);
} else {
return binarySearchSolarLongitude(x, high);
}
}

}
2 changes: 1 addition & 1 deletion base/src/test/java/net/time4j/calendar/hindu/SerializationTest.java
View file
Expand Up @@ -29,7 +29,7 @@ public void serializeVariants() throws IOException, ClassNotFoundException {
roundtrip(rule.variant().with(HinduEra.KALI_YUGA));
roundtrip(rule.variant().withElapsedYears());
roundtrip(rule.variant().withCurrentYears());
roundtrip(rule.variant().withModernHinduSunriseSunset(0.0));
roundtrip(rule.variant().withModernAstronomy(0.0));
roundtrip(rule.variant().withAlternativeLocation(SolarTime.ofMecca()));
}
}
Expand Down

0 comments on commit 451120f

Please sign in to comment.