add a C version

This commit is contained in:
Chen Wei
2014-03-30 16:12:27 +08:00
parent 3e07df2414
commit 73b21aab91
13 changed files with 149952 additions and 0 deletions

26
c/Makefile Normal file
View File

@@ -0,0 +1,26 @@
CC = gcc
CFLAGS = -Wall -O2
LIBS = -lm
SRCS = lunarcal.c astro.c lea406.c vsop.c nutation.c julian.c lunarcalbase.c
OBJS = $(SRCS:.c=.o)
TESTSRCS = testastro.c astro.c lea406.c vsop.c nutation.c julian.c lunarcalbase.c
TESTOBJS = $(TESTSRCS:.c=.o)
LUNARCAL = lunarcal
TESTASTRO = testastro
EXECUTABLE = lunarcal
all: $(LUNARCAL) $(TESTASTRO)
@echo all compiled
$(LUNARCAL): $(OBJS)
$(CC) $(CFLAGS) -o $(LUNARCAL) $(OBJS) $(LIBS)
$(TESTASTRO): $(TESTOBJS)
$(CC) $(CFLAGS) -o $(TESTASTRO) $(TESTOBJS) $(LIBS)
.c.o:
$(CC) $(CFLAGS) -c $< -o $@
clean:
rm -f *.o core a.out astro lunarcal testastro

197
c/astro.c Normal file
View File

@@ -0,0 +1,197 @@
/*
copyright 2014, Chen Wei <weichen302@gmail.com>
version 0.0.3
Implement astronomical algorithms for finding solar terms and moon phases.
Truncated VSOP87D for calculate Sun's apparent longitude;
Truncated LEA-406 for calculate Moon's apparent longitude;
Truncated IAU2000B from USNO NOVAS c source is used for nutation.
Reference:
VSOP87: ftp://ftp.imcce.fr/pub/ephem/planets/vsop87
LEA-406: S. M. Kudryavtsev (2007) "Long-term harmonic development of
lunar ephemeris", Astronomy and Astrophysics 471, 1069-1075
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "astro.h"
#define MAXITER 20 /* max iteration for Secand Method */
/* solve the equation when function f(jd, angle) reaches zero by
* Secand method
*/
double rootbysecand(double (*f)(double , double),
double angle, double x0, double x1, double precision)
{
double fx0, fx1, x2;
fx0 = (*f)(x0, angle);
fx1 = (*f)(x1, angle);
int i = 0;
for (i = 0; i < MAXITER; i++) {
if (fabs(fx1) < precision || fabs(x0 - x1) < precision)
return x1;
x2 = x1 - fx1 * (x1 - x0) / (fx1 - fx0);
fx0 = fx1;
fx1 = (*f)(x2, angle);
x0 = x1;
x1 = x2;
//printf("debug in rootbysecand: iter = %d angle = %.9f \n", i, fx1);
}
return -1;
}
/* covernt radian to 0 - 2pi */
double normrad(double r) {
r = fmod(r, TWOPI);
if (r < 0)
r += TWOPI;
return r;
}
/* convert an angle in radians into (-pi, +pi} */
double npitopi(double r) {
r = fmod(r, TWOPI);
if (r > PI)
r -= TWOPI;
else if (r <= -1.0 * PI)
r += TWOPI;
return r;
}
/* Calculate the difference between target angle and solar geocentric
* longitude at a given JDTT
*
* and normalize the angle between Sun's Longitude on a given
* day and the angle we are looking for to (-pi, pi), therefore f(x) is
* continuous from -pi to pi
*
*/
double f_solarangle(double jd, double angle)
{
return npitopi(apparentsun(jd, 0) - angle);
}
/* Calculate difference between target angle and current sun-moon angle
*
* Arg:
* jd: time in JDTT
* Return:
* angle in radians, convert to -pi to +pi range
*
*/
double f_msangle(double jd, double angle)
{
return npitopi(apparentmoon(jd, 1) - apparentsun(jd, 1) - angle);
}
/* calculate Solar Term by secand method
*
* The Sun's moving speed on ecliptical longitude is 0.04 argsecond / second,
*
* The mean error of truncated VSOP is less than 0.1", nutation by IAU2000B is
* 0.001"
*
* Args:
* year: the year in integer
* angle: degrees of the solar term
* Return:
* time in JDTT
*
*/
double solarterm(int year, double angle)
{
/* mean error when compare apparentsun to NASA(1900-2100) is 0.05"
* 0.000000005 radians = 0.001" */
double ERROR, r, est_vejd, x0, x1;
ERROR = 0.000000005;
GregorianDate ve; /* estimated date of Vernal Equinox, March 20.5 UTC0 */
ve.year = year;
ve.month = 3;
ve.day = 20.5;
est_vejd = g2jd(ve);
/* negative angle means search backward from Vernal Equinox.
* Initialize x0 to the day which apparent Sun longitude close to the
* angle we searching for */
x0 = est_vejd + angle * 360.0 / 365.24;
x1 = x0 + 0.5;
r = angle * DEG2RAD;
return rootbysecand(f_solarangle, r, x0, x1, ERROR);
}
/* search newmoon near a given date.
*
* Angle between Sun-Moon has been converted to {-pi, pi} range so the
* function f_msangle is continuous in that range. Use Secand method to find
* root.
*
* Test shows newmoon can be found in 5 iterations, if the start is close
* enough, it may use only 3 iterations.
*
* Arg:
* jd: in JDTT
* Return:
* JDTT of newmoon
*
*/
double newmoon(double jd)
{
/* 0.0000001 radians is about 0.02 arcsecond, mean error of apparentmoon
* when compared to JPL Horizon is about 0.7 arcsecond */
double ERROR, x0, x1;
ERROR = 0.0000001;
/* initilize x0 to the day close to newmoon */
x0 = jd - f_msangle(jd, 0) / MOON_SPEED;
x1 = x0 + 0.5;
return rootbysecand(f_msangle, 0, x0, x1, ERROR);
}
/* search new moon from specified start time
*
* Arg:
* start: the start time in JD, doesn't matter if it is in TT or UT
* count: the number of newmoons to search after start time
*
* Return:
* a list of JDTT when newmoon occure
*
*/
void findnewmoons(double newmoons[], int nmcount, double startjd)
{
int i;
double nm;
for ( i = 0; i < nmcount; i++) {
nm = newmoon(startjd);
newmoons[i] = nm;
startjd = nm + SYNODIC_MONTH;
}
}
/* convert decimal degree to d m s format string */
size_t fmtdeg(char *strdeg, double d) {
if (abs(d) > 360)
d = fmod(d, 360.0);
double fdegree, deg, m, s, tmp;
fdegree = fabs(d);
tmp = modf(fdegree, &deg);
tmp = modf(tmp * 60.0, &m);
s = tmp * 60;
char *sign = "";
if (d < 0)
sign = "-";
sprintf(strdeg, "%s%d°%d%c%.6f%c", sign, (int)(deg),
(int)(m), 39, s, 34);
return strlen(strdeg);
}

71
c/astro.h Normal file
View File

@@ -0,0 +1,71 @@
/*
* header for astro functions
*/
#define PI 3.14159265358979323846
#define TWOPI 6.28318530717958647693
#define RAD2DEG 57.295779513082320876798
#define DEG2RAD 0.017453292519943295779
#define ASEC2RAD 4.848136811095359935899141e-6
#define ASEC360 1296000.0
#define J2000 2451545.0
#define TROPICAL_YEAR 365.24
#define SYNODIC_MONTH 29.53
#define MOON_SPEED TWOPI / SYNODIC_MONTH /* approximate Moon & Sun's */
#define SUN_SPEED TWOPI / TROPICAL_YEAR /* longitude change per day*/
#define NMCOUNT 15 /* default search total 15 new moons */
#define ISODTLEN 30 /* max length of ISO date string */
typedef struct {
int year;
int month;
double day;
} GregorianDate;
/* Function prototypes */
GregorianDate jd2g(double jd);
size_t fmtdeg(char *strdeg, double d);
double normrad(double r);
double npitopi(double r);
double jdptime(char *isodt, char *fmt, double tz, int isut);
size_t jdftime(char *isodt, double jd, char *fmt, double tz, int isut);
double g2jd(GregorianDate d);
double deltaT(GregorianDate g);
double apparentsun(double jd, int ignorenutation);
double apparentmoon(double jd, int ignorenutation);
double lea406(double jd, int ignorenutation);
double nutation(double jd);
double lightabbr_high(double jd);
double vsopLx(double vsopterms[][3], size_t rowcount, double t);
double vsop(double jd);
double rootbysecand(double (*f)(double , double),
double angle, double x0, double x1, double precision);
double f_solarangle(double jd, double angle);
double f_msangle(double jd, double angle);
double newmoon(double jd);
void findnewmoons(double newmoons[], int nmcount, double startjd);
double solarterm(int year, double angle);
int findastro(int year);

73505
c/jpl_moon.txt Normal file

File diff suppressed because it is too large Load Diff

73505
c/jpl_sun.txt Normal file

File diff suppressed because it is too large Load Diff

368
c/julian.c Normal file
View File

@@ -0,0 +1,368 @@
/*
* functions convert between Gregorian date and Julian date
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "astro.h"
/* convert a Gregorian date to JD from AA, p61 */
double g2jd(GregorianDate d)
{
if (d.month <= 2) {
d.year -= 1;
d.month += 12;
}
int a, b, isjulian;
double jd;
a = (int) (d.year / 100);
isjulian = 0;
if (d.year < 1582) {
isjulian = 1;
} else if (d.year == 1582) {
if (d.month < 10)
isjulian = 1;
if (d.month == 10 && d.day <= 5)
isjulian = 1;
if (d.month == 10 && d.day > 5 && d.day < 15)
return 2299160.5;
}
if (isjulian)
b = 0;
else
b = 2 - a + (int) (a / 4);
/* 30.6001 is a hack Meeus suggested */
jd = (int) (365.25 * (d.year + 4716)) + (int) (30.6001 * (d.month + 1))
+ d.day + b - 1524.5;
return jd;
}
/* convert JD to Gregorian date */
GregorianDate jd2g(double jd)
{
jd += 0.5;
int a, b, c, d, e, z, alpha;
double f;
GregorianDate res;
z = (int) (jd);
f = fmod(jd, 1.0); // decimal part
if (z < 2299161) {
a = z;
} else {
alpha = (int) ((z - 1867216.25) / 36524.25);
a = z + 1 + alpha - (alpha / 4);
}
b = a + 1524;
c = (int) ((b - 122.1) / 365.25);
d = (int) (365.25 * c);
e = (int) ((b - d) / 30.6001);
res.day = (float) (b - d - (int) (30.6001 * e) + f);
if (e < 14)
res.month = e - 1;
else
res.month = e - 13;
if (res.month > 2)
res.year = c - 4716;
else
res.year = c - 4715;
return res;
}
/* Args:
* isodt: datetime string in ISO format
* tz: a integer as timezone, e.g. -8 for UTC-8, 2 for UTC2
* fmt: like strftime but currently only support three format
* %y-%m-%d %H:%M:%S
* %y-%m-%d %H:%M
* %y-%m-%d
* ut: convert to UTC(adjust delta T)
* Return:
* Julian Day
*
*/
double jdptime(char *isodt, char *fmt, double tz, int isut)
{
char inputstr[40];
char *isodate;
char *isot;
strcpy(inputstr, isodt);
if (strcmp(fmt, "%y-%m-%d") == 0) {
char isodatestr[40];
char isotstr[40];
isodate = strcpy(isodatestr, isodt);
isot = strcpy(isotstr, "00:00:00");
} else if (strcmp(fmt, "%y-%m-%d %H:%M") == 0) {
isodate = strtok(inputstr, " ");
isot = strtok(NULL, " ");
} else {
isodate = strtok(inputstr, " ");
isot = strtok(NULL, " ");
}
char *sy;
char *sm, *sd, *shour, *sminute, *ssec;
double d, hour, minute, sec;
GregorianDate g;
sy = strtok(isodate, "-");
sm = strtok(NULL, "-");
sd = strtok(NULL, "-");
g.year = atoi(sy);
g.month = atoi(sm);
d = atof(sd);
shour = strtok(isot, ":");
sminute = strtok(NULL, ":");
ssec = strtok(NULL, ":");
hour = atof(shour);
minute = atof(sminute);
if (ssec)
sec = atof(ssec);
else
sec = 0;
d += (hour * 3600.0 + minute * 60.0 + sec) / 86400.0;
g.day = d;
return g2jd(g);
}
/* format a Julian Day to ISO format datetime
Args:
jd: time in JDTT
tz: a integer as timezone, e.g. -8 for UTC-8, 2 for UTC2
fmt: like strftime but currently only support three format
%y-%m-%d %H:%M:%S
%y-%m-%d %H:%M
%y-%m-%d
ut: convert to UTC(adjust delta T)
Return:
time string in ISO format, e.g. 1984-01-01 23:59:59
*/
size_t jdftime(char *isodt, double jd, char *fmt, double tz, int isut)
{
GregorianDate g;
double deltat, utsec, secs, jdut;
int isecs;
/* char isodt[ISODTLEN]; */
g = jd2g(jd);
deltat = isut ? deltaT(g) : 0;
/* convert jd to seconds, then adjust deltat */
utsec = jd * 86400.0 + tz * 3600.0 - deltat;
jdut = utsec / 86400.0;
/* get time in seconds directly to minimize round error */
secs = fmod(utsec + 43200.0, 86400.0);
if (strcmp(fmt, "%y-%m-%d %H:%M") == 0) {
isecs = (int)(floor(0.5 + secs / 60.0) * 60.0);
} else {
isecs = (int)(secs);
}
if (86400 == secs) {
jdut = floor(jdut) + 0.5;
isecs = 0;
}
g = jd2g(jdut);
/* use integer math hereafter */
int y, m, d, H, M, S, ms;
y = g.year;
m = g.month;
d = floor(g.day);
H = isecs / 3600;
ms = isecs % 3600;
M = ms / 60;
S = ms % 60;
if (strcmp(fmt, "%y-%m-%d") == 0) {
sprintf(isodt, "%04d-%02d-%02d", y, m, d);
} else if (strcmp(fmt, "%y-%m-%d %H:%M") == 0) {
sprintf(isodt, "%04d-%02d-%02d %02d:%02d", y, m, d, H, M);
} else {
sprintf(isodt, "%04d-%02d-%02d %02d:%02d:%02d", y, m, d, H, M, S);
}
return strlen(isodt);
}
/*
Polynomial Expressions for Delta T (ΔT) from nasa. Valid for -1999 to
+3000. http://eclipse.gsfc.nasa.gov/LEcat5/deltatpoly.html
Arg:
year: Gregorian year in integer
m: Gregorian month in integer
d: doesn't matter
Result:
ΔT in seconds
verfify against historical records from NASA
year history computed diff
-500 17190 17195.37 5.4
-400 15530 15523.84 6.2
-300 14080 14071.97 8.0
-200 12790 12786.59 3.4
-100 11640 11632.52 7.5
0 10580 10578.95 1.0
100 9600 9592.45 7.6
200 8640 8636.34 3.7
300 7680 7676.67 3.3
400 6700 6694.67 5.3
500 5710 5705.50 4.5
600 4740 4734.89 5.1
700 3810 3809.04 1.0
800 2960 2951.97 8.0
900 2200 2197.11 2.9
1000 1570 1571.65 1.7
1100 1090 1086.99 3.0
1200 740 735.10 4.9
1300 490 490.98 1.0
1400 320 321.09 1.1
1500 200 197.85 2.2
1600 120 119.55 0.5
1700 9 8.90 0.1
1750 13 13.44 0.4
1800 14 13.57 0.4
1850 7 7.16 0.2
1900 -3 -2.12 0.9
1950 29 29.26 0.3
1955 31 31.23 0.1
1960 33 33.31 0.1
1965 36 36.13 0.4
1970 40 40.66 0.5
1975 46 45.94 0.4
1980 50 50.93 0.4
1985 54 54.60 0.3
1990 57 57.20 0.3
1995 61 61.17 0.4
2000 64 64.00 0.2
2005 65 64.85 0.1
Also, JPL uses "last known leap-second is used over any future interval".
It causes the large error when compare apparent sun/moon position with JPL
Horizon
*/
double deltaT(GregorianDate g) {
int year;
double y, m, u;
year = g.year;
m = (double) g.month;
y = year + (m - 0.5) / 12.0;
if (year < -500) {
u = (year - 1820) / 100.0;
return -20 + 32 * u * u;
} else if (year < 500) {
u = y / 100.0;
return 10583.6 + u * (-1014.41 +
u * (33.78311 +
u * (-5.952053 +
u * (-0.1798452 +
u * (0.022174192 +
u * 0.0090316521)))));
} else if (year < 1600) {
u = (y -1000) / 100.0;
return 1574.2 + u * (-556.01 +
u * (71.23472 +
u * (0.319781 +
u * (-0.8503463 +
u * (-0.005050998 +
u * 0.0083572073)))));
} else if (year < 1700) {
u = y -1600;
return 120 + u * (-0.9808 +
u * (- 0.01532 +
u / 7129));
} else if (year < 1800) {
u = y - 1700;
return 8.83 + u * (0.1603 +
u * (-0.0059285 +
u * (0.00013336 +
u / -1174000)));
} else if (year < 1860) {
u = y - 1800;
return 13.72 + u * (-0.332447 +
u * (0.0068612 +
u * (0.0041116 +
u * (-0.00037436 +
u * (0.0000121272 +
u * (-0.0000001699 +
u * 0.000000000875))))));
} else if (year < 1900) {
u = y - 1860;
return 7.62 + u * (0.5737 +
u * (-0.251754 +
u * (0.01680668 +
u * (-0.0004473624 +
u / 233174))));
} else if (year < 1920) {
u = y - 1900;
return -2.79 + u * (1.494119 +
u * (-0.0598939 +
u * (0.0061966 +
u * -0.000197)));
} else if (year < 1941) {
u = y - 1920;
return 21.20 + u * (0.84493 +
u * (-0.076100 +
u * 0.0020936));
} else if (year < 1961) {
u = y - 1950;
return 29.07 + u * (0.407 +
u * (-1.0 / 233 +
u / 2547));
} else if (year < 1986) {
u = y - 1975;
return 45.45 + u * (1.067 +
u * (-1.0 / 260 +
u / -718));
} else if (year < 2005) {
u = y -2000;
return 63.86 + u * (0.3345 +
u * (-0.060374 +
u * (0.0017275 +
u * (0.000651814 +
u * 0.00002373599))));
/*
* else:
* JPL uses "last known leap-second is used over any future interval" .
* It causes the large error when compare apparent sun/moon position
* with JPL Horizon
*
* return 67.182963
*
*/
} else if (year < 2050) {
u = y -2000;
return 62.92 + u * (0.32217 + u * 0.005589);
} else if (year < 2150) {
u = (y - 1820.0) / 100.0;
return -20 + 32 * u * u - 0.5628 * (2150 - y);
} else {
u = (y - 1820.0) / 100.0;
return -20 + 32 * u * u;
}
}

775
c/lea406.c Normal file
View File

@@ -0,0 +1,775 @@
/*
copyright 2014, Chen Wei <weichen302@gmail.com>
version 0.0.3
Implement astronomical algorithms for finding solar terms and moon phases.
Truncated LEA-406 for calculate Moon's apparent longitude;
Reference:
LEA-406: S. M. Kudryavtsev (2007) "Long-term harmonic development of
lunar ephemeris", Astronomy and Astrophysics 471, 1069-1075
*/
#include <stdio.h>
#include <math.h>
#include "astro.h"
static double FRM[5] = {
785939.924268, 1732564372.3047, -5.279, .006665, -5.522e-5
};
/*
table for LEA-406 moon solution. Those terms are linear combination
of integer multipliers of 14 variables (Arg_j_, j=1,14):
Delaunay variables l, l', F, D;
mean longitude of the ascending node of the Moon {Omega};
mean longitudes of eight major planets {lambda}_pl_;
and the general precession in longitude p_A_.
terms of 3rd-degree and 4th-degree are ignored
in arcsec
average error Moon = 0.73", max 1.5"
*/
static double M_ARG[226][3] = {
{ 485868.249036, 1717915923.21779990, 31.87920 },
{ 1658653.158348, 1488007279.20020032, -44.62040 },
{ 2144521.407384, 3205923202.41800022, -12.74120 },
{ 971736.498072, 3435831846.43559980, 63.75840 },
{ 1287104.793048, 129596581.04809999, -0.55320 },
{ 671559.052464, 3479054525.69560003, -25.50240 },
{ 1172784.909312, -229908644.01759958, -76.49960 },
{ 371548.365300, 1358410698.15210032, -44.06720 },
{ 2630389.656420, 4923839125.63580036, 19.13800 },
{ 857416.614336, 3076326621.36990023, -12.18800 },
{ 801236.544012, -1588319342.16969991, -32.43240 },
{ 1072260.703692, 1602961601.20900011, -6.37060 },
{ 1772973.042084, 1847512504.26589990, 31.32600 },
{ 1472962.354920, -273131323.27759981, 12.76120 },
{ 0.000000, 0.00000000, 0.00000 },
{ 1157427.301500, 5196970448.91339970, 6.37680 },
{ -185690.803428, -1761138602.47780013, 57.38160 },
{ 3803174.565732, 4693930481.61820030, -57.36160 },
{ 1457604.747108, 5153747769.65339947, 95.63760 },
{ 3317306.316696, 2976014558.40040064, -89.24080 },
{ 2945757.951396, 1617603860.24830031, -45.17360 },
{ 3431626.200432, 3335519783.46610022, -13.29440 },
{ 586392.454656, -114954322.00879979, -38.24980 },
{ 2359365.496740, 1732558182.25710011, -6.92380 },
{ 1343284.863372, 4794242544.58769989, 19.69120 },
{ 3116257.905456, 6641755048.85359955, 51.01720 },
{ 2746972.885608, 477658.07790603, -67.34420 },
{ 4289042.814768, 6411846404.83600044, -25.48240 },
{ 686916.660276, -1947824567.23539925, -108.37880 },
{ 315368.294976, -3306235265.38749981, -64.31160 },
{ 2330212.210812, 4967061804.89580059, -70.12280 },
{ -114319.883736, -359505225.06569958, -75.94640 },
{ 1558128.952728, 3320877524.42679977, 25.50860 },
{ -429688.178712, 2946730040.32180023, -11.63480 },
{ 2258841.291120, 3565428427.48369980, 63.20520 },
{ 2574209.586096, 259193162.09619999, -1.10640 },
{ -915556.427748, 1228814117.10400033, -43.51400 },
{ 450160.398036, -6962890.54310000, 7.47222 },
{ 1958830.603956, 1444784599.94020033, 44.64040 },
{ 2816080.459848, 6684977728.11359978, -38.24360 },
{ 2516069.772684, 4564333900.57010078, -56.80840 },
{ 1643295.550536, 6914886372.13119984, 38.25600 },
{ 2730913.862040, 3090968880.40920067, -50.99100 },
{ 3917494.449468, 5053435706.68389988, 18.58480 },
{ 2030201.523648, 2846417977.35230112, -88.68760 },
{ 2088341.337060, -1458722761.12159991, -32.98560 },
{ 4232862.744444, 1747200441.29640031, -45.72680 },
{ 2459889.702360, -100312062.96949959, -77.05280 },
{ 185857.561872, -402727904.32569981, 13.31440 },
{ 4774911.063804, 8129762328.05380058, 6.39680 },
{ 1943472.996144, 6871663692.87119961, 127.51680 },
{ 3001938.021720, 6282249823.78790092, -24.92920 },
{ 100524.205620, -1832870245.22659969, -70.12900 },
{ 2760067.147968, -143534742.22949982, 12.20800 },
{ 300177.445608, -43222679.26000023, 89.26080 },
{ 2845233.745776, 3450474105.47490025, 24.95540 },
{ 2245045.613004, 1373052957.19140100, -82.87020 },
{ 2831438.067660, 1258098635.18260098, -121.12000 },
{ 1829153.112408, 6512158467.80550003, 51.57040 },
{ 3060077.835132, 1977109085.31399989, 30.77280 },
{ 2134680.712596, 1725351443.01107621, -44.04921 },
{ 1873497.247704, 14642259.03930020, -38.80300 },
{ 3602126.154492, 8359670972.07139969, 82.89640 },
{ 3301948.708884, 8402893651.33139992, -6.36440 },
{ 201048.411240, -3665740490.45319939, -140.25800 },
{ 293448.038472, 81066394.15011901, 0.02638 },
{ 1643162.503716, 1717438685.61177826, -210.42628 },
{ 56180.070324, 4664645963.53960037, 20.24440 },
{ 3232841.134644, 1718393581.29570580, -35.46500 },
{ -170499.954060, -5024151188.60529900, -96.19080 },
{ 2538997.273764, 1731206461.65055728, -44.54126 },
{ 238013.777124, 118672081.90543799, 0.28560 },
{ 5090279.358780, 4823527062.66629982, -57.91480 },
{ 2043997.201764, 5038793447.64459991, 57.38780 },
{ 400701.651228, -1876092924.48659992, 19.13180 },
{ 5461827.724080, 6181937760.81840134, -101.98200 },
{ 501225.856848, -3708963169.71319962, -50.99720 },
{ -214844.089356, 1473365020.16090012, -5.81740 },
{ 2744709.540156, 5283344350.70149994, 95.08440 },
{ 1844343.961776, 3249145881.67800045, -102.00200 },
{ 936028.647072, 1710953032.67469978, 39.35142 },
{ -35707.851000, -1724878813.76090002, -24.40698 },
{ -600188.132772, -2077421148.28349924, -107.82560 },
{ 2444698.852992, 3162700523.15799952, 76.51960 },
{ 3963148.815276, 1714428507.27902675, -198.01200 },
{ 1043107.417764, 4837465223.84770012, -69.56960 },
{ 1343118.104928, 6958109051.39120007, -51.00480 },
{ 1958663.845512, 3608651106.74370003, -26.05560 },
{ 3216782.111076, 4808884803.62700081, -19.11180 },
{ 5947695.973116, 7899853684.03620148, -70.10280 },
{ 1528975.666800, 6555381147.06550026, -37.69040 },
{ 1121719.450500, 3472091635.15250015, -18.03018 },
{ 671725.810908, 1315188018.89209986, 45.19360 },
{ 4604411.109744, 3105611139.44850016, -89.79400 },
{ -9960733.695600, 135571.46521378, 63.51699 },
{ 1387628.998668, -1703273664.17849970, -70.68220 },
{ 26423.255340, -32537000.00799108, 44.93644 },
{ -1716792.971760, 2817133459.27370024, -11.08160 },
{ 2129163.799572, 8632802295.34899902, 70.13520 },
{ 1578859.855344, 8210672.20779327, -214.97146 },
{ 1228964.979636, 4434737319.52200031, -56.25520 },
{ 2287389.559392, 162132438.11736214, -200.33773 },
{ 129677.491548, -5067373867.86530018, -6.93000 },
{ 3131615.513268, 1214875955.92260027, -31.85920 },
{ -1401424.676784, -489101806.11379957, -75.39320 },
{ 4975959.475044, 4464021837.60060215, -133.86120 },
{ 4403362.698504, 6771351629.90170002, 50.46400 },
{ 5576147.607816, 6541442985.88409996, -26.03560 },
{ 3487806.270756, 8000165747.00570107, 6.95000 },
{ 4018018.655088, 3220565461.45730019, -51.54420 },
{ 3477280.566240, -3487415.93877248, -229.89120 },
{ 2444532.094548, 5326567029.96150017, 5.82360 },
{ 1743819.756156, 5082016126.90460014, -31.87300 },
{ 2545223.058612, 1329830277.93140078, 6.39060 },
{ 3746994.495408, 29284518.07860041, -77.60600 },
{ 1648812.463560, 7435519.79327631, -75.92841 },
{ -2202661.220796, 1099217536.05590034, -42.96080 },
{ 6176636.963928, 606718.56467720, 146.34632 },
{ 1443809.068992, 2961372299.36110115, -50.43780 },
{ 5260779.312840, 9847678251.27160072, 38.27600 },
{ 1420639.381224, 1369335197.29476237, -44.90600 },
{ 4474733.618196, 8172985007.31380081, -82.86400 },
{ 5299552.000908, 121384894.43375297, 62.45638 },
{ 1602473.088024, -3176638684.33939981, -64.86480 },
{ 476027.554248, 237344163.81087598, 0.57119 },
{ 4718730.993480, 3465116364.51420021, -13.84760 },
{ 9840.694788, 1480571759.40692401, 31.30801 },
{ 987094.105884, -1991047246.49539971, -19.11800 },
{ 1974021.453324, -1818227986.18729925, -108.93200 },
{ -2128477.768848, 113602752.48178008, 40.50488 },
{ 3988865.369160, 6455069084.09600067, -114.74320 },
{ 1906507.630260, 3087251120.51256227, -13.02680 },
{ -4107087.132120, 107747819.74584675, 210.09281 },
{ 863813.842884, 1487527316.94823694, -132.65705 },
{ 628736.401044, -1325874841.08283782, -155.71733 },
{ 602451.478224, -3205445544.34009504, -54.60300 },
{ 2939360.722848, 3206403164.66996431, 75.29545 },
{ 743096.730600, 2716821396.30420065, -88.13440 },
{ -2587922.762544, -1636850170.74401093, 53.56212 },
{ 1714833.228672, 6152653242.73980045, -24.37600 },
{ 2453492.473812, 1488487241.45216417, 43.41625 },
{ 4503886.904124, 4938481384.67510033, -19.66500 },
{ 247854.471912, 1599243841.31236196, 31.59360 },
{ -700712.338392, -244550903.05689979, -37.69660 },
{ -41775.577452, -1650139296.84567857, -64.69155 },
{ 142868.152008, -3043790764.30063820, -187.59653 },
{ 824406.231780, 3722788.71663815, -36.85218 },
{ 3760707.671856, 1406941526.72641158, -130.06172 },
{ -385344.043416, -3550786168.44439936, -102.00820 },
{ -1616186.264472, 1798981675.69158912, 117.32052 },
{ 6433564.222152, 9617769607.25400162, -38.22360 },
{ -394475.866380, 1474716740.76744270, 31.80006 },
{ 3787816.957920, 10120809574.54920006, 25.51480 },
{ 271024.159680, 3191280943.37870026, 26.06180 },
{ -3915532.327356, -1718055041.34457111, -247.79376 },
{ -2943795.829284, 1717776805.09102821, -184.03536 },
{ 4588422.072480, 1333648.70949970, 213.90556 },
{ -327467.049516, -1015133.54049925, 95.02473 },
{ 2429341.245180, 8589579616.08899975, 159.39600 },
{ 723882.026160, 1836588005.12323785, 32.16480 },
{ 3861314.379144, 388789743.14429998, -1.65960 },
{ 1544333.274612, 1128502054.13450146, -120.56680 },
{ 4246575.920892, 3124857449.94421101, -98.18252 },
{ 4168107.377892, 3198959150.64730024, -5.80178 },
{ 1818627.407892, -1491494695.13897252, -185.27080 },
{ 2315021.361444, 8230074391.02330017, 83.44960 },
{ 3331101.994812, 5168390028.69269943, 56.83460 },
{ 15357.607812, -5426879092.93099976, -82.87640 },
{ 10479231.070128, -16418126.94190380, 74.77548 },
{ 1828986.353964, 8676024974.60899925, -19.12560 },
{ 4174722.931032, 6052341179.77030182, -101.42880 },
{ 2577500.752860, -1199777.09773506, -217.20541 },
{ 1182625.604100, 1250663115.38932443, -45.19159 },
{ 123665.467464, 10925660.37799100, -0.30604 },
{ 4960601.867232, 9890900930.53160095, -50.98480 },
{ 3245935.397004, 1574381180.98829985, 44.08720 },
{ 1668493.853136, 2968579038.60712433, -13.31239 },
{ 957940.819956, 1243456376.14330149, -82.31700 },
{ 6461919.320184, 15996023.02879967, 179.55064 },
{ 1101413.989620, -1631542021.42970014, 56.82840 },
{ -377031.566700, -1812274039.13709879, -19.73467 },
{ 1472795.596476, 1890735183.52590013, -57.93480 },
{ -857249.855892, -5240193128.17339993, 82.88400 },
{ 4087994.403528, 10077586895.28919983, 114.77560 },
{ 11070.933228, -1715208268.04166102, 26.29335 },
{ 801403.302456, -3752185848.97319984, 38.26360 },
{ 615545.740584, -3349457944.64750004, 24.94920 },
{ 2014843.915836, 8273297070.28330040, -5.81120 },
{ -2546147.185092, 13289126.10166723, 118.25367 },
{ 4660591.180068, 7770257102.98810196, -69.54960 },
{ 2620548.961632, 3443267366.22887611, -12.17001 },
{ -284819.837796, -5383656413.67099953, -172.13720 },
{ 3008350.635300, -1716566657.86483884, -189.51513 },
{ 3617317.003860, 5096658385.94390011, -70.67600 },
{ 4103352.011340, 4650707802.35820007, 31.89920 },
{ 120935.436876, 3212887254.18870020, -19.68062 },
{ 4103185.252896, 6814574309.16170025, -38.79680 },
{ 3545946.084168, 3695025008.53179979, 62.65200 },
{ 1929677.318028, 4679288222.57890129, -18.55860 },
{ 542048.319360, 6382561886.75740051, 52.12360 },
{ 3682239.128856, 1481043227.42950010, -37.68098 },
{ 1092991.606308, -1709705251.01000595, -246.85066 },
{ 2064728.104380, 1726126595.42559338, -183.09226 },
{ -1101247.231176, -532324485.37379980, 13.86760 },
{ 1310274.480816, 1721633683.11443806, -6.08500 },
{ 1685076.413688, 1455470279.19220901, 0.31604 },
{ -300010.687164, -2120643827.54349971, -18.56480 },
{ -3664055.126520, 2720871.13874409, -227.53801 },
{ 512291.504376, 1685378923.20980859, 76.81564 },
{ 3640898.842560, -1366622384.76644802, 107.07678 },
{ 417145.879260, -1477083922.99626732, -111.06649 },
{ 14236864.720236, -3208579.51268373, -163.76979 },
{ 4389567.020388, 4578976159.60939980, -95.61140 },
{ 786045.694644, 1674693243.95779943, 121.14000 },
{ 150612.259788, 1555781180.92637992, 76.83608 },
{ 459444.993696, 1750452923.22579098, -13.05724 },
{ 176897.182608, 3435351884.18363619, -24.27825 },
{ -3529183.602852, -2591705.71692875, -226.73272 },
{ 169814.944332, -1498930214.93224931, -109.34239 },
{ -356190.757488, -6785289791.08309937, -38.80920 },
{ 3375446.130108, -1329126180.07349992, -33.53880 },
{ 2209276.773936, 1754174550.70710015, -50.44218 },
{ 55466.634672, -1606681665.27969599, -111.04605 },
{ 1131009.331860, -1880047940.86327744, -141.19115 },
{ 6062015.856852, 8259358909.10190010, 5.84360 },
{ 352362.086784, 226418503.43288499, 0.87723 },
};
/*
table for LEA-406 moon solution
Ak0 arcsec Amplitude of the Fourier term
Ak1 mas/yr Amplitude of the 1st-order Poisson term
Ak2 uas/yr2 Amplitude of the 2nd-order Poisson term
{ Ak0, Ak1, Ak2 } ...
*/
static double M_AMP[226][3] = {
{ 22639.5864251, 0.190648, 5.529914 },
{ 4586.4946082, 0.112352, 1.480719 },
{ 2369.9292836, 0.105206, 1.315910 },
{ 769.0251380, 0.013309, 0.373136 },
{ 666.4389699, 16.764899, 0.509567 },
{ 411.5952546, 0.002764, 0.176172 },
{ 211.6566179, 0.006615, 0.018535 },
{ 205.4430995, 5.163882, 0.191171 },
{ 191.9569463, 0.011624, 0.153634 },
{ 164.7319041, 4.138198, 0.167655 },
{ 147.3267606, 3.710258, 0.121729 },
{ 124.9935117, 0.001183, 0.034216 },
{ 109.3839297, 2.755639, 0.083615 },
{ 55.1780671, 0.002891, 0.006919 },
{ 0.0362449, 1.349731, 45.692267 },
{ 45.0995619, 0.000090, 0.030602 },
{ 39.5287366, 0.001225, 0.006971 },
{ 38.4307291, 0.004006, 0.032533 },
{ 36.1241571, 0.000898, 0.025549 },
{ 30.7733692, 0.002140, 0.019268 },
{ 28.3979625, 0.711689, 0.023432 },
{ 24.3591471, 0.612529, 0.022527 },
{ 18.5852985, 0.000437, 0.001521 },
{ 17.9562255, 0.451214, 0.014878 },
{ 14.5306503, 0.364691, 0.018609 },
{ 14.3797216, 0.001391, 0.015092 },
{ 14.2136673, 0.506954, 0.388974 },
{ 13.8992155, 0.002252, 0.015581 },
{ 13.1942141, 0.000494, 0.002232 },
{ 9.6793682, 0.243746, 0.009212 },
{ 9.3659379, 0.000333, 0.007150 },
{ 8.6056849, 0.216461, 0.006770 },
{ 8.4533524, 0.000386, 0.004712 },
{ 8.0506614, 0.404649, 0.009300 },
{ 7.6304201, 0.192240, 0.006663 },
{ 7.4493429, 0.375004, 0.006549 },
{ 7.3717333, 0.370557, 0.007936 },
{ 7.0921386, 0.076778, 0.005079 },
{ 6.3823525, 0.000616, 0.002913 },
{ 5.7416948, 0.000257, 0.005899 },
{ 4.3741497, 0.109740, 0.005590 },
{ 3.9976117, 0.000026, 0.003764 },
{ 3.2097860, 0.000695, 0.001883 },
{ 2.9146498, 0.073322, 0.003241 },
{ 2.7319853, 0.068620, 0.003016 },
{ 2.5682050, 0.129294, 0.002534 },
{ 2.5220218, 0.127230, 0.002163 },
{ 2.4889981, 0.063117, 0.001806 },
{ 2.1461452, 0.053905, 0.001745 },
{ 1.9777716, 0.000440, 0.002642 },
{ 1.9336777, 0.000069, 0.001990 },
{ 1.8708161, 0.046884, 0.002864 },
{ 1.7530780, 0.000160, 0.000246 },
{ 1.4372113, 0.036159, 0.001127 },
{ 1.3713028, 0.000041, 0.000141 },
{ 1.2619335, 0.031743, 0.001156 },
{ 1.2241492, 0.000087, 0.000559 },
{ 1.1868343, 0.000109, 0.000580 },
{ 1.1770227, 0.029544, 0.001701 },
{ 1.1619378, 0.058532, 0.001032 },
{ 1.1431465, 0.001145, 0.004116 },
{ 1.0779824, 0.027072, 0.000890 },
{ 1.0595012, 0.000125, 0.001494 },
{ 0.9902267, 0.000067, 0.001314 },
{ 0.9482810, 0.000049, 0.000381 },
{ 0.8216250, 0.000304, 0.000090 },
{ 0.7791738, 0.055305, 0.066699 },
{ 0.7517241, 0.037784, 0.001034 },
{ 0.7322303, 0.044974, 0.051647 },
{ 0.6694015, 0.016850, 0.000770 },
{ 0.6436261, 0.000332, 0.000422 },
{ 0.6388982, 0.000303, 0.000313 },
{ 0.6352378, 0.016019, 0.000773 },
{ 0.5840142, 0.000053, 0.000496 },
{ 0.5833144, 0.000000, 0.000051 },
{ 0.5715783, 0.000117, 0.000766 },
{ 0.5606404, 0.000000, 0.000176 },
{ 0.5568122, 0.014135, 0.000470 },
{ 0.5459375, 0.013761, 0.000538 },
{ 0.5356650, 0.000028, 0.000298 },
{ 0.5025507, 0.011924, 0.000650 },
{ 0.5005117, 0.011841, 0.000959 },
{ 0.4784073, 0.012033, 0.000371 },
{ 0.4537245, 0.000048, 0.000336 },
{ 0.4381375, 0.054206, 0.022773 },
{ 0.4262242, 0.010710, 0.000506 },
{ 0.4203282, 0.000000, 0.000355 },
{ 0.4134074, 0.010391, 0.000347 },
{ 0.4041956, 0.000444, 0.000404 },
{ 0.3945214, 0.000105, 0.000652 },
{ 0.3821338, 0.009597, 0.000535 },
{ 0.3748102, 0.011708, 0.000597 },
{ 0.3744777, 0.009395, 0.000354 },
{ 0.3575903, 0.009033, 0.000359 },
{ 0.3416863, 0.352069, 0.220915 },
{ 0.3495542, 0.008691, 0.000624 },
{ 0.3435465, 0.009341, 0.007387 },
{ 0.3398415, 0.025631, 0.000333 },
{ 0.3286524, 0.000000, 0.000403 },
{ 0.3248895, 0.001569, 0.031518 },
{ 0.3087352, 0.015508, 0.000464 },
{ 0.3016226, 0.005046, 0.028777 },
{ 0.3015612, 0.007593, 0.000321 },
{ 0.3008704, 0.000080, 0.000190 },
{ 0.2942052, 0.014800, 0.000232 },
{ 0.2925609, 0.000043, 0.000351 },
{ 0.2902265, 0.007303, 0.000385 },
{ 0.2891080, 0.007270, 0.000414 },
{ 0.2825095, 0.007069, 0.000521 },
{ 0.2737910, 0.006879, 0.000263 },
{ 0.2736062, 0.038332, 0.094167 },
{ 0.2633701, 0.006632, 0.000249 },
{ 0.2542937, 0.000000, 0.000192 },
{ 0.2530370, 0.000035, 0.000148 },
{ 0.2500111, 0.012673, 0.000241 },
{ 0.2494213, 0.007693, 0.009143 },
{ 0.2469758, 0.018629, 0.000176 },
{ 0.1180438, 0.242859, 0.232513 },
{ 0.2314117, 0.005794, 0.000255 },
{ 0.2185403, 0.000055, 0.000408 },
{ 0.2111003, 0.000156, 0.000174 },
{ 0.2013358, 0.000026, 0.000300 },
{ 0.1944693, 0.003998, 0.005565 },
{ 0.1931190, 0.009721, 0.000207 },
{ 0.1844919, 0.000278, 0.000754 },
{ 0.1835925, 0.008827, 0.000266 },
{ 0.1825829, 0.000095, 0.000340 },
{ 0.1789891, 0.000029, 0.000060 },
{ 0.1762440, 0.004467, 0.000147 },
{ 0.1751178, 0.006068, 0.003488 },
{ 0.1697756, 0.000016, 0.000207 },
{ 0.1649126, 0.000149, 0.000035 },
{ 0.1643790, 0.002513, 0.016571 },
{ 0.1609152, 0.046251, 0.052193 },
{ 0.1608887, 0.002702, 0.015539 },
{ 0.1603243, 0.024905, 0.034943 },
{ 0.1590996, 0.028613, 0.024866 },
{ 0.1578112, 0.007932, 0.000203 },
{ 0.1537976, 0.004754, 0.006421 },
{ 0.1522670, 0.007645, 0.000263 },
{ 0.1514863, 0.048990, 0.060578 },
{ 0.1499154, 0.003761, 0.000172 },
{ 0.1432907, 0.000164, 0.002489 },
{ 0.1363584, 0.003428, 0.000102 },
{ 0.1360093, 0.005055, 0.007270 },
{ 0.1342991, 0.002202, 0.012985 },
{ 0.1317702, 0.012996, 0.046203 },
{ 0.1312510, 0.004067, 0.005441 },
{ 0.1281174, 0.000034, 0.000093 },
{ 0.1274662, 0.003938, 0.005376 },
{ 0.1261665, 0.000044, 0.000252 },
{ 0.1251903, 0.000032, 0.000034 },
{ 0.1238568, 0.000012, 0.000201 },
{ 0.1207260, 0.003037, 0.000112 },
{ 0.1177733, 0.104057, 0.092375 },
{ 0.1166401, 0.105848, 0.086210 },
{ 0.1145920, 0.065094, 0.062264 },
{ 0.1127967, 0.017578, 0.045202 },
{ 0.1110004, 0.000000, 0.000148 },
{ 0.1100790, 0.000037, 0.000050 },
{ 0.1014638, 0.007699, 0.000024 },
{ 0.0998164, 0.002513, 0.000094 },
{ 0.0992925, 0.003076, 0.004114 },
{ 0.0971982, 0.005531, 0.000173 },
{ 0.0935175, 0.010642, 0.005973 },
{ 0.0932035, 0.002337, 0.000160 },
{ 0.0920486, 0.002316, 0.000087 },
{ 0.0915362, 0.000000, 0.000057 },
{ 0.0909223, 0.014292, 0.002950 },
{ 0.0909178, 0.000000, 0.000101 },
{ 0.0903302, 0.002261, 0.000156 },
{ 0.0585008, 0.089697, 0.078271 },
{ 0.0891047, 0.000146, 0.000132 },
{ 0.0860207, 0.003064, 0.003710 },
{ 0.0849998, 0.000015, 0.000149 },
{ 0.0847363, 0.002131, 0.000072 },
{ 0.0840293, 0.000147, 0.000081 },
{ 0.0831051, 0.002088, 0.000051 },
{ 0.0829235, 0.005841, 0.006803 },
{ 0.0828709, 0.002083, 0.000063 },
{ 0.0824031, 0.004719, 0.002548 },
{ 0.0805004, 0.002028, 0.000058 },
{ 0.0801760, 0.000019, 0.000000 },
{ 0.0776552, 0.000012, 0.000133 },
{ 0.0755097, 0.005767, 0.017169 },
{ 0.0752280, 0.000083, 0.000049 },
{ 0.0750004, 0.001874, 0.000066 },
{ 0.0737291, 0.001851, 0.000119 },
{ 0.0718922, 0.005260, 0.006839 },
{ 0.0714181, 0.001784, 0.000140 },
{ 0.0686150, 0.000124, 0.000210 },
{ 0.0684983, 0.000000, 0.000045 },
{ 0.0676977, 0.013106, 0.022804 },
{ 0.0674175, 0.001690, 0.000069 },
{ 0.0659954, 0.000011, 0.000077 },
{ 0.0657981, 0.003704, 0.000148 },
{ 0.0654081, 0.001645, 0.000080 },
{ 0.0651462, 0.003283, 0.000055 },
{ 0.0650640, 0.001653, 0.000084 },
{ 0.0643910, 0.003236, 0.000104 },
{ 0.0643077, 0.003626, 0.000118 },
{ 0.0640864, 0.000767, 0.007024 },
{ 0.0633090, 0.000680, 0.006279 },
{ 0.0631287, 0.003177, 0.000038 },
{ 0.0622945, 0.001764, 0.000987 },
{ 0.0617621, 0.001669, 0.001381 },
{ 0.0610645, 0.001547, 0.000186 },
{ 0.0610286, 0.034491, 0.021496 },
{ 0.0605601, 0.001685, 0.001334 },
{ 0.0590447, 0.001206, 0.001683 },
{ 0.0576070, 0.003983, 0.003844 },
{ 0.0574169, 0.012520, 0.037732 },
{ 0.0572562, 0.000022, 0.000070 },
{ 0.0567906, 0.000000, 0.000000 },
{ 0.0561596, 0.005382, 0.001305 },
{ 0.0559131, 0.001479, 0.001084 },
{ 0.0530823, 0.005285, 0.000393 },
{ 0.0527616, 0.030032, 0.013528 },
{ 0.0525526, 0.004791, 0.003973 },
{ 0.0516483, 0.001301, 0.000065 },
{ 0.0514262, 0.003899, 0.000000 },
{ 0.0513864, 0.003162, 0.000101 },
{ 0.0510076, 0.005738, 0.003400 },
{ 0.0507210, 0.001707, 0.002544 },
{ 0.0507028, 0.001274, 0.000084 },
{ 0.0506960, 0.005853, 0.006825 },
};
/*
table for LEA-406 moon solution
phik0 arcsec Phase of the Fourier term
phik1 arcsec Phase of the 1st-order Poisson term
phik2 arcsec Phase of the 2nd-order Poisson term
{ phik0, phik1, phik2 } ...
*/
static double M_PHASE[226][3] = {
{ 0.000383901368, -85.864631587374, -90.051320519700 },
{ 0.002465002781, 43.853451435402, -89.760808801732 },
{ 0.002701282720, 9.198419860616, -89.927193714594 },
{ 0.000733167975, -69.133923809631, -89.854596360650 },
{ 179.998822844683, 0.033165685600, -4.896329847950 },
{ 179.999821873246, -92.231132704615, 90.039281835354 },
{ 0.002182942141, 115.589973490862, -89.724570062688 },
{ 0.004492444068, 179.946598681046, -152.700235600240 },
{ 0.003074613545, 3.237177076528, -89.045873461368 },
{ 0.004755315788, 179.925373896526, -139.458634581667 },
{ 179.998402135869, 0.029599610744, -21.201860289276 },
{ -179.998700951455, -21.042848422574, 85.696683781596 },
{ 179.998830243184, 0.038692795619, 12.231960211961 },
{ 0.002824379245, -1.618145480441, -93.712463437426 },
{ 90.000000000000, -90.000000000000, -90.000000000000 },
{ -179.999730775095, 136.948333197864, 90.252913650594 },
{ 0.000631501778, -27.748331128161, 91.027469629070 },
{ 0.004768684748, 15.338642008082, -91.448391559281 },
{ 0.000827624832, -44.890207530000, -94.464735632788 },
{ 0.004679600181, 28.129778606743, -91.253712037737 },
{ 179.997410623704, -0.012385396825, 19.711675596544 },
{ -179.999633530725, -0.005264403413, 33.706044239720 },
{ -179.998546923256, -52.237586501500, 94.104800351413 },
{ -0.001202833150, -179.943767466285, -161.488197277814 },
{ 0.006060440686, 179.946114884688, -128.325040292843 },
{ 0.003479681319, -5.088467767454, -88.879853308440 },
{ 76.222761080242, 131.248177879490, 158.937351156777 },
{ 0.005491991867, 5.559956085120, -89.285634473759 },
{ 0.001670679988, 119.052640751395, 95.344742980084 },
{ 179.997512146320, 0.012469626353, -38.622239081360 },
{ -179.997591940179, -130.753370044026, 90.132938283259 },
{ 0.004027243751, 179.894863506410, -167.685307438556 },
{ -179.997888374255, -164.568228146555, 94.313669562461 },
{ 0.004280309127, 179.942703717388, -140.168177475912 },
{ 179.999403180215, 0.057104471446, 29.127407493447 },
{ -179.986595625958, -0.009533448592, -7.291531705498 },
{ 0.004516315551, 179.946856948035, -147.990210285414 },
{ -2.292723477875, 179.558678822974, 99.917394355623 },
{ -179.995813706528, 174.359278613627, 85.899138385237 },
{ -179.997028119590, -170.875555577117, 87.517830497342 },
{ 0.007586565777, 179.860744031295, -126.955336135073 },
{ -179.999173757590, 173.529013740000, 89.645385356640 },
{ -179.996333121479, -173.529824098564, 90.691511437192 },
{ -179.998440539648, 0.000498095926, 46.506659248843 },
{ 0.007351480495, 179.846498962627, -134.320841941915 },
{ 179.998955202477, 0.004715888323, -24.817267057930 },
{ -179.974500785427, 0.037906703343, 16.516954102688 },
{ 0.004084722791, 179.949205415543, -177.339668887924 },
{ 0.005169265951, 179.947259672081, -162.782325495969 },
{ 0.005470571565, 7.802235215572, -91.188350231592 },
{ 0.001979936864, -49.093247644398, -89.853610040000 },
{ 0.008430819626, 179.862846789852, -120.579155305642 },
{ -179.997429267858, 3.689592432343, -8.984125112071 },
{ 179.999682180708, 0.034556477923, 0.880098853934 },
{ -179.998734104791, 23.889956013343, 92.315478473154 },
{ -0.001089631169, -179.953391970037, -148.465693395845 },
{ -179.995711702960, -44.744675254758, 91.037547248403 },
{ 0.005618410521, 141.891652358833, -89.245762243791 },
{ 0.006176245341, 179.923324075525, -121.671975846595 },
{ -179.990710441378, -0.006065350786, 4.820672054727 },
{ -179.850243125450, 10.554417037974, 75.954784760505 },
{ 0.004538058638, 179.991538999959, -178.668525637314 },
{ 0.004948942002, -2.305853377327, -89.055962190853 },
{ -179.996220310084, -162.841071063425, 90.115626432671 },
{ 0.001362391462, 124.153219100000, 89.375195054863 },
{ -179.987972335590, -9.455484053495, 98.486655980000 },
{ 100.520850457896, 58.245822583006, 176.656756595668 },
{ 0.005558909986, 179.939402800171, -129.032487586412 },
{ 74.866340599937, 113.143692095032, 157.140609294154 },
{ 179.996802948398, 0.023597286958, -49.820634167458 },
{ 179.991295527846, 105.627702205222, 104.722979510000 },
{ 1.229379639312, 14.658638231892, 174.161916553362 },
{ -179.997302559910, -0.079892424376, 50.903582526549 },
{ -179.997054101519, -179.406547096258, 88.325519523174 },
{ -179.997923631717, -21.306014230000, -111.306014230000 },
{ 0.008920446563, 16.507710861469, -87.001504751198 },
{ -179.997496718333, -49.229404680000, -96.386960454652 },
{ 179.993029626649, 0.063953608450, 32.857661723889 },
{ 179.999568497340, 0.029889128486, 41.670863171014 },
{ -179.997678361057, -131.554420507361, 90.182772169547 },
{ -6.998845124057, -179.836968150987, 95.494728149239 },
{ -6.946998653071, -179.864579320978, 76.439431277611 },
{ 0.003950080511, 179.877986481689, 178.009646557383 },
{ -179.995214967737, -178.270591376070, 92.820667680820 },
{ -169.915008504888, 163.606608570999, 141.706947652864 },
{ -179.995168304093, -0.129249467750, 51.249058012792 },
{ -0.000588147800, 76.911637520000, -89.818659700943 },
{ -0.001541311372, -179.960228699023, -157.839192649548 },
{ 0.005446540202, 178.120205971300, -91.077675513499 },
{ 0.009961559799, 12.474727024070, -91.103784962745 },
{ -179.994627410236, -0.125960849016, 58.152530864523 },
{ -10.452602212154, -175.909739396951, 95.972435843391 },
{ -179.993551894629, -0.075640934119, 40.075819425085 },
{ -179.999730936646, -0.111110307792, 44.048546642547 },
{ -77.540885610051, -152.586445555125, 101.407916729860 },
{ 0.012860009915, 179.181823683989, -38.176426666159 },
{ -159.862611524727, 8.333796508938, 108.553241211148 },
{ 0.000894638748, 179.956445219452, -108.623909313499 },
{ -179.998497415727, -141.434388770000, 89.925582532280 },
{ 33.792050127818, 148.766561177540, 127.241668563970 },
{ 0.009118370918, 179.900634928826, -123.775932072104 },
{ -112.134619005092, -20.559070831297, -20.121839387134 },
{ -0.002033212332, -179.939721200975, 132.145526586787 },
{ 0.006967986668, 11.254792408294, -89.887013647003 },
{ 0.004247032092, -179.958621060283, -154.368735202973 },
{ 0.009442611117, 35.703533951771, -89.937681650476 },
{ -179.997723052131, -0.017116494452, 54.758385850839 },
{ -179.995947079637, -0.068525925581, 57.217045528267 },
{ 0.009819658662, 179.838703270416, -112.860829392226 },
{ 0.001459400774, 179.958199617562, -139.757079434325 },
{ -167.517828096945, 168.948506773399, 132.172642029344 },
{ -0.001169845673, -179.992214129080, -139.036692930000 },
{ 0.001621751548, -34.394376710000, -89.706775576900 },
{ -179.994253678680, -163.461615332369, 90.508788209221 },
{ 179.954707818255, 0.297463641368, 79.938876374848 },
{ -0.288803957182, 138.401774201477, -90.194124627437 },
{ -0.000441652594, 179.983728062620, -113.490522358660 },
{ 113.159919693549, 21.205321636711, -68.950861756151 },
{ -179.993471286898, -0.138018683023, 48.941925280000 },
{ 0.008513287422, 9.383107446051, -89.982476550838 },
{ 178.266197023425, 167.040019546567, 63.337092639167 },
{ -179.993411147698, -152.143374643636, 90.067296178865 },
{ 177.709033274638, -94.278194341766, 86.467689692292 },
{ 179.998212230083, 0.011725483834, -38.211204526568 },
{ -179.548460175794, 11.748806769316, 155.794165646379 },
{ 179.696182246530, 8.388992907636, -155.810091409325 },
{ -179.966541471463, 123.169858230588, 174.061712674745 },
{ 0.001277362909, -87.646654344739, 28.795532201940 },
{ 0.006441557128, 179.856257333904, 171.660707410000 },
{ -73.684186429799, 66.397350353519, -166.301305668267 },
{ -179.993990076579, -136.968231614654, 90.111372676863 },
{ 178.195523468679, 167.215437886350, 100.414547150000 },
{ 33.641453360113, -165.302705060093, -60.601265278652 },
{ -69.150402372136, -59.249946857851, -144.807244596642 },
{ -112.036976032111, -23.538862057265, -19.703539149818 },
{ 71.225252170846, 89.953326699449, 149.522242168418 },
{ -106.906299285279, -128.752837957579, -35.167482742815 },
{ 0.008709768323, 179.926987734485, -128.528021306430 },
{ 125.434998080106, -144.834917590238, 33.658347760855 },
{ 0.009988580138, 179.896870091933, -116.342563520000 },
{ -112.167984339780, -118.315481678004, -31.324502874434 },
{ 0.001914516745, 179.958829517227, -128.259077640388 },
{ 178.919893127688, 111.151535580000, 127.496580553665 },
{ 179.997049090149, 0.013601633772, -2.620698856245 },
{ -72.167265771598, -161.189710332737, 19.930112286947 },
{ -111.985346564208, -23.928553663926, -20.132412817928 },
{ 29.213939211730, 160.784938812168, -84.562976989473 },
{ 54.598102587697, -34.751396609636, 146.388967571845 },
{ 179.997627140643, -20.956099913023, -107.488668694704 },
{ 125.438781294500, -144.945781138477, 33.725620650288 },
{ 0.011052840008, 12.898827180000, -89.856260770306 },
{ 179.874812111593, -138.253779548140, 27.781164827339 },
{ -179.995364958510, -152.171377200000, 90.034521688113 },
{ 179.999571716845, 0.088736680056, 43.415194343448 },
{ 15.485979928550, 79.333070815266, 171.905198795860 },
{ 15.814692146435, 81.813245835498, 173.563281819956 },
{ 131.151069984933, -154.395177439830, -38.414622359212 },
{ -44.655851925818, 117.757239601126, 162.851093530505 },
{ 0.002900038346, 135.182987450000, -89.948686156302 },
{ 1.103187588210, 22.976396007451, -111.078340600000 },
{ -179.887622646628, -0.016270483931, -31.472745223488 },
{ 0.007797056081, 179.864439387698, -138.165985101854 },
{ 54.595537656712, -34.659546623063, 146.299996591100 },
{ -80.604321141847, 23.349211093147, 124.690639100328 },
{ 9.950565354210, -20.646007710451, -20.232712520916 },
{ 0.007228561268, 179.913696450803, -114.674710042095 },
{ 0.000031181430, -179.969483846305, -115.306109670000 },
{ -179.998385048108, 85.733997830000, -94.266002170000 },
{ 126.785821299905, 46.152057795628, 40.582093829691 },
{ 0.000067308052, -58.051764990000, -90.033259843394 },
{ 0.011879132634, 179.748408306343, -114.695467227172 },
{ -170.230856074017, 100.037571146737, -1.882374403492 },
{ -0.844779089926, 176.042515812187, -118.079402338350 },
{ -73.380351208713, 159.697723721737, -39.602387509092 },
{ -179.992838936462, -160.881735064112, 90.089950323242 },
{ 0.008160895360, 179.745911260715, -149.577199957246 },
{ -0.946094322429, 178.710437293623, -77.509104380747 },
{ -179.993591655201, -0.214491674788, 93.905327790000 },
{ 156.746173873619, -64.861968042177, 49.372967047162 },
{ -179.999469648778, -0.016408343945, -35.948330450000 },
{ 37.928299901718, 129.004508305459, 132.215335260474 },
{ 179.998321758011, 0.020516890633, 4.395108701959 },
{ -179.991344977137, 26.519098943619, -31.875040030000 },
{ 0.005727141458, 4.355481318395, -90.149163090774 },
{ -21.479998361998, -146.605834790916, 58.779750145468 },
{ -179.964083052415, -55.875067383234, 47.387971540000 },
{ -0.001625399610, -179.904773285909, 142.248731360116 },
{ -179.993620429096, -0.123104591073, 70.321134490000 },
{ -162.315246851023, -72.633139513377, 112.510164638260 },
{ 0.012390392562, 179.760835926114, -111.438742365499 },
{ -179.872748658637, 50.117769415320, 168.302109917160 },
{ 0.000466648467, 169.116621610000, 79.116621610000 },
{ 68.455744968383, -114.610536287778, -179.414411109053 },
{ -0.002100097173, 179.947053742654, -135.363706761591 },
{ -179.991954083944, -149.820003150000, 89.975569316466 },
{ -102.019559358387, 161.184331716547, 85.320435106082 },
{ 0.000547819855, 179.949477302242, -126.016097727692 },
{ -179.987915461393, -0.051829095371, 5.014976620001 },
{ 0.008347581988, 179.849649173293, -123.589509971061 },
{ 0.006800628510, 179.936140204095, -119.801617824531 },
{ -73.109646466056, 9.949188683703, 119.696875238960 },
{ 38.312689959023, -174.350323749453, 129.275979873543 },
{ 38.530542970871, -178.317385228219, 132.123003310646 },
{ 0.006989175447, -179.665485193957, -124.587860108589 },
{ 177.277845106886, 126.415631259050, -109.793598391998 },
{ -161.442918729888, 2.240961856056, 106.127634621147 },
{ -179.990868095964, -3.697578545908, 7.251882193224 },
{ 61.035930743232, -40.663545438817, -108.996557846408 },
{ -158.629498383529, 16.851256552374, 110.695920982574 },
{ 176.824703041161, -98.747691569943, 85.257139856799 },
{ 177.683627303808, -120.178114678752, -84.678375180739 },
{ 145.813265306575, 165.236834102475, -115.686032935051 },
{ -179.990848592847, -174.767338145272, 87.971105506103 },
{ 179.998967305008, -128.346026290000, -128.346026290000 },
{ 109.706873518703, -161.692320544537, 8.071930386460 },
{ -21.373532585359, -175.453622486338, 97.727650370240 },
{ -76.964196141885, -25.052062779649, -88.184334546093 },
{ 174.044963900841, 69.697068381257, -103.438549174491 },
{ 1.411566523810, -63.763692231235, 96.605146342888 },
{ -0.002484367922, -179.922758430555, 124.122955805312 },
{ -179.975293067009, 0.144068585072, -127.623925030000 },
{ -85.595401423341, 33.309416345973, 52.614419166008 },
{ -2.682922479904, 90.346879348731, 101.825176323036 },
{ -72.194222741362, -164.226260049548, 24.011727797009 },
{ -179.994957071833, -0.078650599201, 64.442567813218 },
{ -167.705956663408, -121.678544311070, 135.639537931150 },
};
/*
* LEA-406 Moon Solution
*
* Reference:
* Long-term harmonic development of lunar ephemeris.
* Kudryavtsev S.M. <Astron. Astrophys. 471, 1069 (2007)>
*
*/
/* compute moon ecliptic longitude using lea406 */
double lea406(double jd, int ignorenutation) {
double t, tm, tm2;
t = (jd - J2000) / 36525.0;
tm = t / 10.0;
tm2 = tm * tm;
double V, arg;
V = FRM[0] + (((FRM[4] * t + FRM[3]) * t + FRM[2]) * t + FRM[1]) * t;
int i;
for (i = 0; i < 226; i++) {
arg = (M_ARG[i][0] + t * (M_ARG[i][1] + M_ARG[i][2] * t)) * ASEC2RAD;
V += M_AMP[i][0] * sin(arg + M_PHASE[i][0] * DEG2RAD)
+ M_AMP[i][1] * sin(arg + M_PHASE[i][1] * DEG2RAD) * tm
+ M_AMP[i][2] * sin(arg + M_PHASE[i][2] * DEG2RAD) * tm2;
}
V *= ASEC2RAD;
if (!ignorenutation) {
V += nutation(jd);
/* printf("debug lea406, nutation been adjusted"); */
}
return V;
}
/* calculate the apparent position of the Moon, it is an alias to the
* lea406 function
*/
double apparentmoon(double jd, int ignorenutation)
{
return lea406(jd, ignorenutation);
}

26
c/lunarcal.c Normal file
View File

@@ -0,0 +1,26 @@
#include <stdio.h>
#include <stdlib.h>
#include "lunarcalbase.h"
int main(int argc, char *argv[])
{
int start, end;
if (argc == 2) {
start = atoi(argv[1]);
end = start;
} else if (argc == 3) {
start = atoi(argv[1]);
end = atoi(argv[2]);
} else {
printf("Usage: lunarcal startyear endyear \n");
exit(2);
}
while (start <= end) {
cn_lunarcal(start);
start++;
}
return 0;
}

276
c/lunarcalbase.c Normal file
View File

@@ -0,0 +1,276 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "astro.h"
#include "lunarcalbase.h"
static char *CN_DAY[] = {
"", "",
"初二", "初三", "初四", "初五", "初六", "初七", "初八", "初九",
"初十", "十一", "十二", "十三", "十四", "十五", "十六", "十七",
"十八", "十九", "二十", "廿一", "廿二", "廿三", "廿四", "廿五",
"廿六", "廿七", "廿八", "廿九", "三十"
};
static char *CN_MON[] = {
"十二月",
"正月", "二月", "三月", "四月", "五月", "六月",
"七月", "八月", "九月", "十月", "十一月", "十二月", /* index = 12 */
"", "", "", "", "", "", /* pad from index 13 to 18 */
"閏十一月", "閏十二月", "閏正月", "閏二月", "閏三月", "閏四月",
"閏五月", "閏六月", "閏七月", "閏八月", "閏九月", "閏十月",
"閏十一月", "閏十二月"
};
/* solar term's angle + 120 / 15 */
static char *CN_SOLARTERM[] = {
"小雪", "大雪", "冬至", "小寒", "大寒", "立春", "雨水", "驚蟄",
"春分", "清明", "穀雨", "立夏", "小滿", "芒種", "夏至", "小暑",
"大暑", "立秋", "處暑", "白露", "秋分", "寒露", "霜降", "立冬",
"小雪", "大雪", "冬至"
};
static double newmoons[MAX_NEWMOONS];
static double solarterms[MAX_SOLARTERMS];
static struct lunarcal *lcs[MAX_DAYS];
static int firstnm_offset;
static int cached_year[CACHESIZE]; /* caches intermedia lunar cal*/
static struct lunarcal *cached_lcs[CACHESIZE][MAX_DAYS];
int cachep = 0; /* next free location in cache */
/* normalize Julian Day to midnight after adjust timezone and deltaT */
double normjd(double jd, double tz)
{
double deltat, tmp, jd_i, jd_f;
GregorianDate g;
g = jd2g(jd);
deltat = deltaT(g);
tmp = jd + (tz * 3600.0 - deltat) / 86400.0;
g = jd2g(tmp);
jd_f = modf(tmp, &jd_i);
if (jd_f >= 0.5)
jd = jd_i + 0.5;
else
jd = jd_i - 0.5;
return jd;
}
void cn_lunarcal(int year)
{
int i, k, len1, len2;
double ystart, yend;
struct lunarcal *thisyear[MAX_DAYS];
struct lunarcal *nextyear[MAX_DAYS];
struct lunarcal *output[MAX_DAYS];
len1 = get_cached_lc(thisyear, year);
len2 = get_cached_lc(nextyear, year + 1);
/* zip */
GregorianDate g;
g.year = year;
g.month = 1;
g.day = 1;
ystart = g2jd(g);
g.month = 12;
g.day = 31;
yend = g2jd(g);
k = 0;
for (i = 0; i < MAX_DAYS; output[i++] = NULL)
;
for (i = 0; i < len1; i++) {
if (thisyear[i]->jd >= ystart)
output[k++] = thisyear[i];
}
for (i = 0; i < len2; i++) {
if (nextyear[i]->jd <= yend)
output[k++] = nextyear[i];
}
print_lunarcal(output, k);
}
int get_cache_index(int year)
{
int i;
for (i = 0; i < CACHESIZE; i++) {
if (cached_year[i] == year)
return i;
}
return -1;
}
int get_cached_lc(struct lunarcal *p[], int year)
{
int i, k, len;
double angle;
double jd_nm, est_nm;
if ((k = get_cache_index(year)) > -1) {
for (i = 0; i < MAX_DAYS; i++) {
if (cached_lcs[k][i] ==NULL)
break;
p[i] = cached_lcs[k][i];
}
return i;
}
for (i = 0; i < MAX_SOLARTERMS; i++) {
angle = (double) (i * 15 - 120);
solarterms[i] = normjd(solarterm(year, angle), 8);
}
/* search 15 newmoons start 30 days before last Winter Solstice */
est_nm = solarterms[2] - 30;
for (i = 0; i < MAX_NEWMOONS; i++) {
jd_nm = newmoon(est_nm);
newmoons[i] = normjd(jd_nm, 8);
est_nm = jd_nm + SYNODIC_MONTH;
}
len = mark_month_day();
for (i = 0; i < MAX_DAYS; i++) {
if (cachep > CACHESIZE)
cachep = 0;
cached_lcs[cachep][i] = lcs[i];
p[i] = lcs[i];
}
/* add to cache */
cached_year[cachep++] = year;
return len;
}
/* mark month and day number, solarterms */
int mark_month_day(void)
{
int i, k, len;
int leapmonth, month;
double lc_start, lc_end, jd, month_day1;
struct lunarcal *lc;
for (k = 0; k < MAX_DAYS; lcs[k++] = NULL)
;
month = 11;
leapmonth = find_leap();
lc_start = newmoons[firstnm_offset];
lc_end = solarterms[MAX_SOLARTERMS -1];
jd = lc_start;
len = 0;
while (jd < lc_end) {
/* scan for month jd belongs */
for (i = firstnm_offset; i < MAX_NEWMOONS; i++) {
if (jd < newmoons[i]) {
month = 11 + i - 1 - firstnm_offset;
month_day1 = newmoons[i - 1];
break;
}
}
/* adjust leapmonth */
if (leapmonth > -1.0 && month == leapmonth) {
month = (month - 1) % 12 + 20; /* add offset 20 to distinguish */
} else if (leapmonth > -1.0 && month > leapmonth) {
month = (month - 1) % 12;
} else {
month %= 12;
}
lc = lcalloc(jd);
lc->month = month;
lc->day = (int) (jd - month_day1 + 1);
lcs[len++] = lc;
jd += 1.0;
}
/* add solar terms to tree */
for (i = 0; i < MAX_SOLARTERMS - 1; i++) {
if (solarterms[i] >= lc_start) {
lc = lcalloc(solarterms[i]);
k = (int) (lc->jd - lc_start);
if (lcs[k])
lcs[k]->solarterm = i;
}
}
return len;
}
int find_leap(void)
{
/* count newmoons between two Winter Solstice */
int nmcount, flag, leapmonth, i, n;
nmcount = 0;
for (i = 0; i < MAX_NEWMOONS; i++) {
if (newmoons[i] > solarterms[2]
&& newmoons[i] <= solarterms[MAX_SOLARTERMS - 1])
nmcount += 1;
}
/* first new moon of new year */
for (i = 0; i < MAX_NEWMOONS; i++) {
if (newmoons[i] > solarterms[2]) {
firstnm_offset = i - 1;
break;
}
}
/* leap year has more than 12 newmoons between two Winter Solstice */
leapmonth = -1;
if (nmcount > 12) {
for (i = 1; i < MAX_NEWMOONS; i++) {
flag = 0;
for (n = 0; n < MAX_SOLARTERMS; n++) {
if (solarterms[n] >= newmoons[i - 1]
&& solarterms[n] < newmoons[i]
&& n % 2 == 0)
flag = 1;
}
if (!flag) {
leapmonth = 11 + i - 1 - firstnm_offset;
break;
}
}
}
return leapmonth;
}
struct lunarcal *lcalloc(double jd)
{
struct lunarcal *p;
p = (struct lunarcal *) malloc(sizeof(struct lunarcal));
if (p) {
p->jd = jd;
p->solarterm = -1;
p->month = -1;
p->day = -1;
}
return p;
}
void print_lunarcal(struct lunarcal *p[], int len)
{
int i;
char isodate[30];
char tmp[30];
for (i = 0; i < len; i++) {
jdftime(isodate, p[i]->jd, "%y-%m-%d", 0, 0);
if (p[i]->day == 1)
sprintf(tmp, "%s", CN_MON[p[i]->month]);
else
sprintf(tmp, "%s", CN_DAY[p[i]->day]);
if (p[i]->solarterm > -1)
sprintf(tmp, "%s %s", tmp, CN_SOLARTERM[p[i]->solarterm]);
printf("%s %s \n", isodate, tmp);
}
}

30
c/lunarcalbase.h Normal file
View File

@@ -0,0 +1,30 @@
#define MAX_SOLARTERMS 27
#define MAX_NEWMOONS 15
#define MAX_DAYS 450
#define CACHESIZE 5
struct lunarcal {
double jd;
int solarterm; /* index of solarterm */
int month; /* month in Lunar Calendar */
int day; /* day in Lunar Calendar */
};
/* Function prototypes */
void cn_lunarcal(int year);
int get_cached_lc(struct lunarcal *p[], int year);
double normjd(double jd, double tz);
int find_leap(void);
int mark_month_day(void);
struct lunarcal *lcalloc(double jd);
void print_lunarcal(struct lunarcal *p[], int len);
int get_cache_index(int year);

220
c/nutation.c Normal file
View File

@@ -0,0 +1,220 @@
/*
copyright 2014, Chen Wei <weichen302@gmail.com>
version 0.0.3
Implement astronomical algorithms for finding solar terms and moon phases.
Truncated IAU2000B from USNO NOVAS c source is used for nutation.
The higher accuracy light abberation algorithm from A & A. p156,
*/
#include <stdio.h>
#include <math.h>
#include "astro.h"
/*
77 terms IAU2000B Luni-Solar nutation table from USNO NOVAS c source
Luni-Solar argument multipliers, coefficients, unit 1e-7 arcsec
L L' F D Om longitude (sin, t*sin, cos), obliquity (cos, t*cos, sin)
*/
static double IAU2000BNutationTable[77][11] = {
{ 0, 0, 0, 0, 1, -172064161, -174666, 33386, 92052331, 9086, 15377 },
{ 0, 0, 2, -2, 2, -13170906, -1675, -13696, 5730336, -3015, -4587 },
{ 0, 0, 2, 0, 2, -2276413, -234, 2796, 978459, -485, 1374 },
{ 0, 0, 0, 0, 2, 2074554, 207, -698, -897492, 470, -291 },
{ 0, 1, 0, 0, 0, 1475877, -3633, 11817, 73871, -184, -1924 },
{ 0, 1, 2, -2, 2, -516821, 1226, -524, 224386, -677, -174 },
{ 1, 0, 0, 0, 0, 711159, 73, -872, -6750, 0, 358 },
{ 0, 0, 2, 0, 1, -387298, -367, 380, 200728, 18, 318 },
{ 1, 0, 2, 0, 2, -301461, -36, 816, 129025, -63, 367 },
{ 0, -1, 2, -2, 2, 215829, -494, 111, -95929, 299, 132 },
{ 0, 0, 2, -2, 1, 128227, 137, 181, -68982, -9, 39 },
{ -1, 0, 2, 0, 2, 123457, 11, 19, -53311, 32, -4 },
{ -1, 0, 0, 2, 0, 156994, 10, -168, -1235, 0, 82 },
{ 1, 0, 0, 0, 1, 63110, 63, 27, -33228, 0, -9 },
{ -1, 0, 0, 0, 1, -57976, -63, -189, 31429, 0, -75 },
{ -1, 0, 2, 2, 2, -59641, -11, 149, 25543, -11, 66 },
{ 1, 0, 2, 0, 1, -51613, -42, 129, 26366, 0, 78 },
{ -2, 0, 2, 0, 1, 45893, 50, 31, -24236, -10, 20 },
{ 0, 0, 0, 2, 0, 63384, 11, -150, -1220, 0, 29 },
{ 0, 0, 2, 2, 2, -38571, -1, 158, 16452, -11, 68 },
{ 0, -2, 2, -2, 2, 32481, 0, 0, -13870, 0, 0 },
{ -2, 0, 0, 2, 0, -47722, 0, -18, 477, 0, -25 },
{ 2, 0, 2, 0, 2, -31046, -1, 131, 13238, -11, 59 },
{ 1, 0, 2, -2, 2, 28593, 0, -1, -12338, 10, -3 },
{ -1, 0, 2, 0, 1, 20441, 21, 10, -10758, 0, -3 },
{ 2, 0, 0, 0, 0, 29243, 0, -74, -609, 0, 13 },
{ 0, 0, 2, 0, 0, 25887, 0, -66, -550, 0, 11 },
{ 0, 1, 0, 0, 1, -14053, -25, 79, 8551, -2, -45 },
{ -1, 0, 0, 2, 1, 15164, 10, 11, -8001, 0, -1 },
{ 0, 2, 2, -2, 2, -15794, 72, -16, 6850, -42, -5 },
{ 0, 0, -2, 2, 0, 21783, 0, 13, -167, 0, 13 },
{ 1, 0, 0, -2, 1, -12873, -10, -37, 6953, 0, -14 },
{ 0, -1, 0, 0, 1, -12654, 11, 63, 6415, 0, 26 },
{ -1, 0, 2, 2, 1, -10204, 0, 25, 5222, 0, 15 },
{ 0, 2, 0, 0, 0, 16707, -85, -10, 168, -1, 10 },
{ 1, 0, 2, 2, 2, -7691, 0, 44, 3268, 0, 19 },
{ -2, 0, 2, 0, 0, -11024, 0, -14, 104, 0, 2 },
{ 0, 1, 2, 0, 2, 7566, -21, -11, -3250, 0, -5 },
{ 0, 0, 2, 2, 1, -6637, -11, 25, 3353, 0, 14 },
{ 0, -1, 2, 0, 2, -7141, 21, 8, 3070, 0, 4 },
{ 0, 0, 0, 2, 1, -6302, -11, 2, 3272, 0, 4 },
{ 1, 0, 2, -2, 1, 5800, 10, 2, -3045, 0, -1 },
{ 2, 0, 2, -2, 2, 6443, 0, -7, -2768, 0, -4 },
{ -2, 0, 0, 2, 1, -5774, -11, -15, 3041, 0, -5 },
{ 2, 0, 2, 0, 1, -5350, 0, 21, 2695, 0, 12 },
{ 0, -1, 2, -2, 1, -4752, -11, -3, 2719, 0, -3 },
{ 0, 0, 0, -2, 1, -4940, -11, -21, 2720, 0, -9 },
{ -1, -1, 0, 2, 0, 7350, 0, -8, -51, 0, 4 },
{ 2, 0, 0, -2, 1, 4065, 0, 6, -2206, 0, 1 },
{ 1, 0, 0, 2, 0, 6579, 0, -24, -199, 0, 2 },
{ 0, 1, 2, -2, 1, 3579, 0, 5, -1900, 0, 1 },
{ 1, -1, 0, 0, 0, 4725, 0, -6, -41, 0, 3 },
{ -2, 0, 2, 0, 2, -3075, 0, -2, 1313, 0, -1 },
{ 3, 0, 2, 0, 2, -2904, 0, 15, 1233, 0, 7 },
{ 0, -1, 0, 2, 0, 4348, 0, -10, -81, 0, 2 },
{ 1, -1, 2, 0, 2, -2878, 0, 8, 1232, 0, 4 },
{ 0, 0, 0, 1, 0, -4230, 0, 5, -20, 0, -2 },
{ -1, -1, 2, 2, 2, -2819, 0, 7, 1207, 0, 3 },
{ -1, 0, 2, 0, 0, -4056, 0, 5, 40, 0, -2 },
{ 0, -1, 2, 2, 2, -2647, 0, 11, 1129, 0, 5 },
{ -2, 0, 0, 0, 1, -2294, 0, -10, 1266, 0, -4 },
{ 1, 1, 2, 0, 2, 2481, 0, -7, -1062, 0, -3 },
{ 2, 0, 0, 0, 1, 2179, 0, -2, -1129, 0, -2 },
{ -1, 1, 0, 1, 0, 3276, 0, 1, -9, 0, 0 },
{ 1, 1, 0, 0, 0, -3389, 0, 5, 35, 0, -2 },
{ 1, 0, 2, 0, 0, 3339, 0, -13, -107, 0, 1 },
{ -1, 0, 2, -2, 1, -1987, 0, -6, 1073, 0, -2 },
{ 1, 0, 0, 0, 2, -1981, 0, 0, 854, 0, 0 },
{ -1, 0, 0, 1, 0, 4026, 0, -353, -553, 0, -139 },
{ 0, 0, 2, 1, 2, 1660, 0, -5, -710, 0, -2 },
{ -1, 0, 2, 4, 2, -1521, 0, 9, 647, 0, 4 },
{ -1, 1, 0, 1, 1, 1314, 0, 0, -700, 0, 0 },
{ 0, -2, 2, -2, 1, -1283, 0, 0, 672, 0, 0 },
{ 1, 0, 2, 2, 1, -1331, 0, 8, 663, 0, 4 },
{ -2, 0, 2, 2, 2, 1383, 0, -2, -594, 0, -2 },
{ -1, 0, 0, 0, 2, 1405, 0, 4, -610, 0, 2 },
{ 1, 1, 2, -2, 2, 1290, 0, 0, -556, 0, 0 }
};
/*
Calculate nutation angles using the IAU 2000B model.
The table is from NOVAS(USNO) c source flle.
"...the 77-term, IAU-approved truncated nutation series, IAU 2000B, which
is accurate to about 0.001 arcsecond in the interval 1995-2050"
Arg:
jd as jd
Return:
nutation of longitude in radians
*/
double nutation(double jd) {
double t, L, Lp, F, D, Om;
t = (jd - J2000) / 36525.0;
/* Mean anomaly of the Moon, in arcsec */
L = 485868.249036 + t * 1717915923.2178;
/* Mean anomaly of the Sun. */
Lp = 1287104.79305 + t * 129596581.0481;
/* Mean argument of the latitude of the Moon. */
F = 335779.526232 + t * 1739527262.8478;
/* Mean elongation of the Moon from the Sun. */
D = 1072260.70369 + t * 1602961601.2090;
/* Mean longitude of the ascending node of the Moon. */
Om = 450160.398036 - t * 6962890.5431;
double arg, lon, deplan;
lon = 0;
int i;
for (i = 0; i < 77; i++) {
arg = ( IAU2000BNutationTable[i][0] * L
+ IAU2000BNutationTable[i][1] * Lp
+ IAU2000BNutationTable[i][2] * F
+ IAU2000BNutationTable[i][3] * D
+ IAU2000BNutationTable[i][4] * Om) * ASEC2RAD;
lon += ( IAU2000BNutationTable[i][5]
+ IAU2000BNutationTable[i][6] * t) * sin(arg)
+ IAU2000BNutationTable[i][7] * cos(arg);
}
/* unit of longitude is 1.0e-7 arcsec, convert it to arcsec */
lon *= 1.0e-7;
/* Constant account for the missing long-period planetary terms in the
truncated nutation model, in arcsec */
deplan = 0.000388;
lon += deplan;
lon *= ASEC2RAD;
/* lon = fmod(lon, TWOPI); */
return lon;
}
/* the higher accuracy light abberation algorithm from A & A. p156,
* the error will be less than 0.001"
*/
double lightabbr_high(double jd)
{
double tm, tm2, tm3, t, t2, t3;
t = (jd - J2000) / 36525.0;
t2 = t * t;
t3 = t * t2;
tm = t / 10.0;
tm2 = tm * tm;
tm3 = tm * tm2;
// variation of the Sun's longitude
double var_lon;
var_lon = 3548.330
+ 118.568 * sin(1.527664004990 + 6283.075850600876 * tm)
+ 2.476 * sin(1.484508993906 + 12566.151699456424 * tm)
+ 1.376 * sin(0.486077687339 + 77713.771468687730 * tm)
+ 0.119 * sin(1.276490181677 + 7860.419392621536 * tm)
+ 0.114 * sin(5.885711004647 + 5753.384884566103 * tm)
+ 0.086 * sin(3.884055717388 + 11506.769769132206 * tm)
+ 0.078 * sin(2.841633387025 + 161000.685738008644 * tm)
+ 0.054 * sin(1.441333038870 + 18849.227550057301 * tm)
+ 0.052 * sin(2.993569534399 + 3930.209696310768 * tm)
+ 0.034 * sin(0.529208263814 + 71430.695618086858 * tm)
+ 0.033 * sin(2.091087703461 + 5884.926847413107 * tm)
+ 0.023 * sin(4.320419446313 + 5223.693920276658 * tm)
+ 0.023 * sin(5.674983441420 + 5507.553238331428 * tm)
+ 0.021 * sin(2.707426294193 + 11790.629088932305 * tm)
+ 7.311 * tm * sin(5.819826570714 + 6283.075850600876 * tm)
+ 0.305 * tm * sin(5.776715192860 + 12566.151699456424 * tm)
+ 0.010 * tm * sin(5.733703298774 + 18849.227550057301 * tm)
+ 0.309 * tm2 * sin(4.214128894867 + 6283.075850600876 * tm)
+ 0.021 * tm2 * sin(3.578766215288 + 12566.151699456424 * tm)
+ 0.004 * tm3 * sin(5.198655163283 + 77713.771468687730 * tm)
+ 0.010 * tm3 * sin(2.700139544566 + 6283.075850600876 * tm);
double M, e, C, v, R, res;
// mean anomaly of the Sun
M = (357.52910 + 35999.0503 * t - 0.0001559 * t2
- 0.00000048 * t3) * DEG2RAD;
// the eccentricity of Earth's orbit
e = 0.016708617 - 0.000042037 * t - 0.0000001236 * t2;
// Sun's equation of center
C = ( (1.9146 - 0.004817 * t - 0.000014 * t2) * sin(M)
+ (0.019993 - 0.000101 * t) * sin(2 * M)
+ 0.00029 * sin(3 * M)) * DEG2RAD;
// true anomaly
v = M + C;
// Sun's distance from the Earth, in AU
R = (1.000001018 * (1 - e * e)) / (1 + e * cos(v));
res = -0.005775518 * R * var_lon * ASEC2RAD;
return res;
}

237
c/testastro.c Normal file
View File

@@ -0,0 +1,237 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "astro.h"
#define MAX_JPL_LINE_LEN 100
#define MAX_JPL_RECORDS 73415
#define FLAG_SOE 1
struct jplrcd {
double jd;
double lon;
};
int parsejplhorizon(char *fname, struct jplrcd *records[]);
struct jplrcd *lon_alloc(void);
void verify_apparent_sun_moon(void);
double n180to180(double angle);
double jd2year(double jd);
double jd2year(double jd)
{
double fyear, jdyearstart;
GregorianDate tmp, yearstart ;
tmp = jd2g(jd);
yearstart.year = tmp.year;
yearstart.month = 1;
yearstart.day = 1;
jdyearstart = g2jd(yearstart);
fyear = (double) tmp.year + (jd - jdyearstart) / 365.0;
return fyear;
}
void testdeltat(void);
void testdeltat()
{
// double d = -133.5;
// int i;
double jd;
char strout[30];
jd = jdptime("2012-01-05 18:00", "%y-%m-%d %H:%M", 0, 0);
//jd = jdptime("2012-01-05", "%y-%m-%d", 0, 0);
printf("%f\n", jd);
jdftime(strout, jd, "%y-%m-%d %H:%M", 0, 0);
printf("jdftime output = %s\n", strout);
GregorianDate g;
int i,year;
double deltat;
year = -500;
for (i = 0; i < 20; i++) {
g.year = year;
g.month = 1;
g.day = 0;
deltat = deltaT(g);
printf("%d = %.2f\n", g.year, deltat);
year += 100;
}
return;
}
void testnewmoon_solarterm(void);
void testnewmoon_solarterm(void)
{
double newmoons[NMCOUNT];
double jd;
//jd = jdptime("2014-01-01 18:00", "%y-%m-%d %H:%M", 0, 0);
GregorianDate g;
int year = 2000;
int n;
char isodt[30];
int i;
for (n = 0; n < 5; n++) {
g.year = year;
g.month = 1;
g.day = 1;
jd = g2jd(g);
findnewmoons(newmoons, NMCOUNT, jd);
year += 1;
}
for (i = 0; i < NMCOUNT; i++) {
jdftime(isodt, newmoons[i], "%y-%m-%d %H:%M:%S", 8.0, 1);
printf("found newmoon: %s %.8f\n", isodt, newmoons[i]);
}
double angle;
for (angle = -90; angle < 285; angle += 15) {
jd = solarterm(2014, angle);
jdftime(isodt, jd, "%y-%m-%d %H:%M:%S", 8.0, 1);
printf("solar term: %3.0f %s\n", angle, isodt);
}
return;
}
void testapparentmoon(void);
void testapparentmoon(void)
{
double jd = 2411545.0;
char deg[30];
char degsun[30];
double d;
int i;
for (i = 0; i < 20; i++) {
d = apparentmoon(jd, 1) * RAD2DEG;
fmtdeg(deg, d);
d = lightabbr_high(jd) * RAD2DEG;
fmtdeg(degsun, d);
printf("%.2f %s %s\n", jd, deg, degsun);
jd += 2;
}
}
void testnutation(void);
void testnutation(void)
{
double jd = 2411545.0;
char deg[30];
double d;
d = nutation(jd) * RAD2DEG;
fmtdeg(deg, d);
printf("%s\n", deg);
}
static char buf[MAX_JPL_LINE_LEN];
int parsejplhorizon(char *fname, struct jplrcd *records[])
{
FILE *fp;
char *p;
int i, flag;
struct jplrcd *plon;
if ((fp = fopen(fname, "r")) == NULL) {
printf("can not open %s\n", fname);
exit(2);
}
flag = 0;
i = 0;
while ((p = fgets(buf, MAX_JPL_LINE_LEN + 1, fp)) != NULL
&& i < MAX_JPL_RECORDS) {
if (p == strstr(p, "$$SOE")) { /* start of records */
flag = 1;
continue;
} else if (p == strstr(p, "$$EOE")) {
break;
}
if (flag) {
plon = lon_alloc();
sscanf(p, "%lf %lf", &(plon->jd), &(plon->lon));
records[i++] = plon;
}
}
return i;
}
struct jplrcd *lon_alloc(void)
{
return (struct jplrcd *) malloc(sizeof(struct jplrcd));
}
/* verify accuracy against JPL
* output can be save and used for gnuplot
*/
void verify_apparent_sun_moon(void)
{
int lensun, lenmoon;
int i, step, count;
double delta_sun, delta_moon;
double delta_sun_n, delta_sun_p, delta_moon_n, delta_moon_p;
struct jplrcd *jplsun[MAX_JPL_RECORDS];
struct jplrcd *jplmoon[MAX_JPL_RECORDS];
lensun = parsejplhorizon("jpl_sun.txt", jplsun);
lenmoon = parsejplhorizon("jpl_moon.txt", jplmoon);
step = 50;
i = 0;
count = 0;
delta_sun_n = 0;
delta_sun_p = 0;
delta_moon_n = 0;
delta_moon_p = 0;
while (i < lensun) {
if (jplsun[i]->jd == jplmoon[i]->jd) {
delta_sun = n180to180(apparentsun(jplsun[i]->jd, 0) * RAD2DEG
- jplsun[i]->lon) * 3600;
delta_moon = n180to180(apparentmoon(jplmoon[i]->jd, 0) * RAD2DEG
- jplmoon[i]->lon) * 3600;
if (delta_sun > 0)
delta_sun_p += delta_sun;
else
delta_sun_n += delta_sun;
if (delta_moon > 0)
delta_moon_p += delta_moon;
else
delta_moon_n += delta_moon;
count++;
printf("%.2f %.9f %.9f\n",
jd2year(jplsun[i]->jd), delta_moon, delta_sun);
}
i += step;
}
printf("\n# total records of JPL Sun = %d Moon=%d\n", lensun, lenmoon);
printf("# Mean Error (arcsec):\n");
printf("# Sun: +%.4f / %.4f Moon: +%.4f / %.4f\n",
delta_sun_p / count, delta_sun_n / count,
delta_moon_p / count, delta_moon_n / count);
}
double n180to180(double angle)
{
angle = fmod(angle, 360.0);
if (angle > 180.0)
angle += -360.0;
else if (angle <= -180.0)
angle += 360;
return angle;
}
int main(void);
int main()
{
//testnewmoon_solarterm();
//testapparentmoon();
//testnutation();
verify_apparent_sun_moon();
return 0;
}

716
c/vsop.c Normal file
View File

@@ -0,0 +1,716 @@
/*
copyright 2014, Chen Wei <weichen302@gmail.com>
version 0.0.3
Implement astronomical algorithms for finding solar terms and moon phases.
Truncated VSOP87D for calculate Sun's apparent longitude;
Reference:
VSOP87: ftp://ftp.imcce.fr/pub/ephem/planets/vsop87
*/
#include <stdio.h>
#include <math.h>
#include "astro.h"
/* Truncated VSOP87D tables */
static double earth_L0[172][3] = {
{ 1.75347045673, 0, 0 },
{ 0.03341656456, 4.66925680417, 6283.0758499914 },
{ 0.00034894275, 4.62610241759, 12566.1516999828 },
{ 0.00003417571, 2.82886579606, 3.523118349 },
{ 0.00003497056, 2.74411800971, 5753.3848848968 },
{ 0.00003135896, 3.62767041758, 77713.7714681205 },
{ 0.00002676218, 4.41808351397, 7860.4193924392 },
{ 0.00002342687, 6.13516237631, 3930.2096962196 },
{ 0.00001273166, 2.03709655772, 529.6909650946 },
{ 0.00001324292, 0.74246356352, 11506.7697697936 },
{ 0.00000901855, 2.04505443513, 26.2983197998 },
{ 0.00001199167, 1.10962944315, 1577.3435424478 },
{ 0.00000857223, 3.50849156957, 398.1490034082 },
{ 0.00000779786, 1.17882652114, 5223.6939198022 },
{ 0.0000099025, 5.23268129594, 5884.9268465832 },
{ 0.00000753141, 2.53339053818, 5507.5532386674 },
{ 0.00000505264, 4.58292563052, 18849.2275499742 },
{ 0.00000492379, 4.20506639861, 775.522611324 },
{ 0.00000356655, 2.91954116867, 0.0673103028 },
{ 0.00000284125, 1.89869034186, 796.2980068164 },
{ 0.0000024281, 0.34481140906, 5486.777843175 },
{ 0.00000317087, 5.84901952218, 11790.6290886588 },
{ 0.00000271039, 0.31488607649, 10977.078804699 },
{ 0.0000020616, 4.80646606059, 2544.3144198834 },
{ 0.00000205385, 1.86947813692, 5573.1428014331 },
{ 0.00000202261, 2.45767795458, 6069.7767545534 },
{ 0.00000126184, 1.0830263021, 20.7753954924 },
{ 0.00000155516, 0.83306073807, 213.299095438 },
{ 0.00000115132, 0.64544911683, 0.9803210682 },
{ 0.00000102851, 0.63599846727, 4694.0029547076 },
{ 0.00000101724, 4.26679821365, 7.1135470008 },
{ 0.00000099206, 6.20992940258, 2146.1654164752 },
{ 0.00000132212, 3.41118275555, 2942.4634232916 },
{ 0.00000097607, 0.6810127227, 155.4203994342 },
{ 0.00000085128, 1.29870743025, 6275.9623029906 },
{ 0.00000074651, 1.75508916159, 5088.6288397668 },
{ 0.00000101895, 0.97569221824, 15720.8387848784 },
{ 0.00000084711, 3.67080093025, 71430.69561812909 },
{ 0.00000073547, 4.67926565481, 801.8209311238 },
{ 0.00000073874, 3.50319443167, 3154.6870848956 },
{ 0.00000078756, 3.03698313141, 12036.4607348882 },
{ 0.00000079637, 1.807913307, 17260.1546546904 },
{ 0.00000085803, 5.98322631256, 161000.6857376741 },
{ 0.00000056963, 2.78430398043, 6286.5989683404 },
{ 0.00000061148, 1.81839811024, 7084.8967811152 },
{ 0.00000069627, 0.83297596966, 9437.762934887 },
{ 0.00000056116, 4.38694880779, 14143.4952424306 },
{ 0.00000062449, 3.97763880587, 8827.3902698748 },
{ 0.00000051145, 0.28306864501, 5856.4776591154 },
{ 0.00000055577, 3.47006009062, 6279.5527316424 },
{ 0.00000041036, 5.36817351402, 8429.2412664666 },
{ 0.00000051605, 1.33282746983, 1748.016413067 },
{ 0.00000051992, 0.18914945834, 12139.5535091068 },
{ 0.00000049, 0.48735065033, 1194.4470102246 },
{ 0.000000392, 6.16832995016, 10447.3878396044 },
{ 0.00000035566, 1.77597314691, 6812.766815086 },
{ 0.0000003677, 6.04133859347, 10213.285546211 },
{ 0.00000036596, 2.56955238628, 1059.3819301892 },
{ 0.00000033291, 0.59309499459, 17789.845619785 },
{ 0.00000035954, 1.70876111898, 2352.8661537718 },
{ 0.00000040938, 2.39850881707, 19651.048481098 },
{ 0.00000030047, 2.73975123935, 1349.8674096588 },
{ 0.00000030412, 0.44294464135, 83996.84731811189 },
{ 0.00000023663, 0.48473567763, 8031.0922630584 },
{ 0.00000023574, 2.06527720049, 3340.6124266998 },
{ 0.00000021089, 4.14825464101, 951.7184062506 },
{ 0.00000024738, 0.21484762138, 3.5904286518 },
{ 0.00000025352, 3.16470953405, 4690.4798363586 },
{ 0.0000002282, 5.22197888032, 4705.7323075436 },
{ 0.00000021419, 1.42563735525, 16730.4636895958 },
{ 0.00000021891, 5.55594302562, 553.5694028424 },
{ 0.00000017481, 4.56052900359, 135.0650800354 },
{ 0.00000019925, 5.22208471269, 12168.0026965746 },
{ 0.0000001986, 5.77470167653, 6309.3741697912 },
{ 0.000000203, 0.37133792946, 283.8593188652 },
{ 0.00000014421, 4.19315332546, 242.728603974 },
{ 0.00000016225, 5.98837722564, 11769.8536931664 },
{ 0.00000015077, 4.19567181073, 6256.7775301916 },
{ 0.00000019124, 3.82219996949, 23581.2581773176 },
{ 0.00000018888, 5.38626880969, 149854.40013480789 },
{ 0.00000014346, 3.72355084422, 38.0276726358 },
{ 0.00000017898, 2.21490735647, 13367.9726311066 },
{ 0.00000012054, 2.62229588349, 955.5997416086 },
{ 0.00000011287, 0.17739328092, 4164.311989613 },
{ 0.00000013971, 4.40138139996, 6681.2248533996 },
{ 0.00000013621, 1.88934471407, 7632.9432596502 },
{ 0.00000012503, 1.13052412208, 5.5229243074 },
{ 0.00000010498, 5.35909518669, 1592.5960136328 },
{ 0.00000009803, 0.99947478995, 11371.7046897582 },
{ 0.0000000922, 4.57138609781, 4292.3308329504 },
{ 0.00000010327, 6.19982566125, 6438.4962494256 },
{ 0.00000012003, 1.003514567, 632.7837393132 },
{ 0.00000010827, 0.32734520222, 103.0927742186 },
{ 0.00000008356, 4.53902685948, 25132.3033999656 },
{ 0.00000010005, 6.0291496328, 5746.271337896 },
{ 0.00000008409, 3.29946744189, 7234.794256242 },
{ 0.00000008006, 5.82145271907, 28.4491874678 },
{ 0.00000010523, 0.93871805506, 11926.2544136688 },
{ 0.00000007686, 3.12142363172, 7238.6755916 },
{ 0.00000009378, 2.62414241032, 5760.4984318976 },
{ 0.00000008127, 6.11228001785, 4732.0306273434 },
{ 0.00000009232, 0.48343968736, 522.5774180938 },
{ 0.00000009802, 5.24413991147, 27511.4678735372 },
{ 0.00000007871, 0.99590177926, 5643.1785636774 },
{ 0.00000008123, 6.2705301365, 426.598190876 },
{ 0.00000009048, 5.33686335897, 6386.16862421 },
{ 0.0000000862, 4.16538210888, 7058.5984613154 },
{ 0.00000006297, 4.71724819317, 6836.6452528338 },
{ 0.00000007575, 3.97382858911, 11499.6562227928 },
{ 0.00000007756, 2.95729056763, 23013.5395395872 },
{ 0.00000007314, 0.60652505806, 11513.8833167944 },
{ 0.00000005955, 2.87641047971, 6283.14316029419 },
{ 0.00000006534, 5.79072926033, 18073.7049386502 },
{ 0.00000007188, 3.99831508699, 74.7815985673 },
{ 0.00000007346, 4.38582365437, 316.3918696566 },
{ 0.00000005413, 5.39199024641, 419.4846438752 },
{ 0.00000005127, 2.36062848786, 10973.55568635 },
{ 0.00000007056, 0.32258441903, 263.0839233728 },
{ 0.00000006625, 3.66475158672, 17298.1823273262 },
{ 0.00000006762, 5.91132535899, 90955.5516944961 },
{ 0.00000004938, 5.73672165674, 9917.6968745098 },
{ 0.00000005547, 2.45152597661, 12352.8526045448 },
{ 0.00000005958, 3.32051344676, 6283.0085396886 },
{ 0.00000004471, 2.06385999536, 7079.3738568078 },
{ 0.00000006153, 1.45823331144, 233141.31440436149 },
{ 0.00000004348, 4.4234217548, 5216.5803728014 },
{ 0.00000006123, 1.07494905258, 19804.8272915828 },
{ 0.00000004488, 3.6528503715, 206.1855484372 },
{ 0.0000000402, 0.83995823171, 20.3553193988 },
{ 0.00000005188, 4.06503864016, 6208.2942514241 },
{ 0.00000005307, 0.38217636096, 31441.6775697568 },
{ 0.00000003785, 2.34369213733, 3.881335358 },
{ 0.00000004497, 3.27230796845, 11015.1064773348 },
{ 0.00000004132, 0.92128915753, 3738.761430108 },
{ 0.00000003521, 5.97844807108, 3894.1818295422 },
{ 0.00000004215, 1.90601120623, 245.8316462294 },
{ 0.00000003701, 5.03069397926, 536.8045120954 },
{ 0.00000003865, 1.82634360607, 11856.2186514245 },
{ 0.00000003652, 1.01838584934, 16200.7727245012 },
{ 0.0000000339, 0.97785123922, 8635.9420037632 },
{ 0.00000003737, 2.95380107829, 3128.3887650958 },
{ 0.00000003507, 3.71291946325, 6290.1893969922 },
{ 0.00000003086, 3.64646921512, 10.6366653498 },
{ 0.00000003397, 1.10590684017, 14712.317116458 },
{ 0.00000003334, 0.83684924911, 6496.3749454294 },
{ 0.00000002805, 2.58504514144, 14314.1681130498 },
{ 0.0000000365, 1.08344142571, 88860.05707098669 },
{ 0.00000003388, 3.20185096055, 5120.6011455836 },
{ 0.00000003252, 3.47859752062, 6133.5126528568 },
{ 0.00000002553, 3.94869034189, 1990.745017041 },
{ 0.0000000352, 2.05559692878, 244287.60000722769 },
{ 0.00000002565, 1.560717849, 23543.23050468179 },
{ 0.00000002621, 3.85639359951, 266.6070417218 },
{ 0.00000002955, 3.39692949667, 9225.539273283 },
{ 0.00000002876, 6.02635617464, 154717.60988768269 },
{ 0.00000002395, 1.16131956403, 10984.1923516998 },
{ 0.00000003161, 1.32798718453, 10873.9860304804 },
{ 0.00000003163, 5.08946464629, 21228.3920235458 },
{ 0.00000002361, 4.27212906992, 6040.3472460174 },
{ 0.0000000303, 1.80209931347, 35371.8872659764 },
{ 0.00000002343, 3.576898605, 10969.9652576982 },
{ 0.00000002618, 2.57870156528, 22483.84857449259 },
{ 0.00000002113, 3.71393780256, 65147.6197681377 },
{ 0.00000002019, 0.81393923319, 170.6728706192 },
{ 0.00000002003, 0.38091017375, 6172.869528772 },
{ 0.00000002506, 3.74379142438, 10575.4066829418 },
{ 0.00000002381, 0.10581361289, 7.046236698 },
{ 0.00000001949, 4.86892513469, 36.0278666774 },
{ 0.00000002074, 4.2279477457, 5650.2921106782 },
{ 0.00000001924, 5.5946054986, 6282.0955289232 },
{ 0.00000001949, 1.07002512703, 5230.807466803 },
{ 0.00000001988, 5.19736046771, 6262.300454499 },
};
static double earth_L1[165][3] = {
{ 6283.31966747491, 0, 0 },
{ 0.00206058863, 2.67823455584, 6283.0758499914 },
{ 0.0000430343, 2.63512650414, 12566.1516999828 },
{ 0.00000425264, 1.59046980729, 3.523118349 },
{ 0.00000108977, 2.96618001993, 1577.3435424478 },
{ 0.00000093478, 2.59212835365, 18849.2275499742 },
{ 0.00000119261, 5.79557487799, 26.2983197998 },
{ 0.00000072122, 1.13846158196, 529.6909650946 },
{ 0.00000067768, 1.87472304791, 398.1490034082 },
{ 0.00000067327, 4.40918235168, 5507.5532386674 },
{ 0.00000059027, 2.8879703846, 5223.6939198022 },
{ 0.00000055976, 2.17471680261, 155.4203994342 },
{ 0.00000045407, 0.39803079805, 796.2980068164 },
{ 0.00000036369, 0.46624739835, 775.522611324 },
{ 0.00000028958, 2.64707383882, 7.1135470008 },
{ 0.00000019097, 1.84628332577, 5486.777843175 },
{ 0.00000020844, 5.34138275149, 0.9803210682 },
{ 0.00000018508, 4.96855124577, 213.299095438 },
{ 0.00000016233, 0.03216483047, 2544.3144198834 },
{ 0.00000017293, 2.99116864949, 6275.9623029906 },
{ 0.00000015832, 1.43049285325, 2146.1654164752 },
{ 0.00000014615, 1.20532366323, 10977.078804699 },
{ 0.00000011877, 3.25804815607, 5088.6288397668 },
{ 0.00000011514, 2.07502418155, 4694.0029547076 },
{ 0.00000009721, 4.23925472239, 1349.8674096588 },
{ 0.00000009969, 1.30262991097, 6286.5989683404 },
{ 0.00000009452, 2.69957062864, 242.728603974 },
{ 0.00000012461, 2.83432285512, 1748.016413067 },
{ 0.00000011808, 5.2737979048, 1194.4470102246 },
{ 0.00000008577, 5.64475868067, 951.7184062506 },
{ 0.00000010641, 0.76614199202, 553.5694028424 },
{ 0.00000007576, 5.30062664886, 2352.8661537718 },
{ 0.00000005834, 1.76649917904, 1059.3819301892 },
{ 0.00000006385, 2.65033984967, 9437.762934887 },
{ 0.00000005223, 5.66135767624, 71430.69561812909 },
{ 0.00000005305, 0.90857521574, 3154.6870848956 },
{ 0.00000006101, 4.66632584188, 4690.4798363586 },
{ 0.0000000433, 0.24102555403, 6812.766815086 },
{ 0.00000005041, 1.42490103709, 6438.4962494256 },
{ 0.00000004259, 0.77355900599, 10447.3878396044 },
{ 0.00000005198, 1.85353197345, 801.8209311238 },
{ 0.00000003744, 2.00119516488, 8031.0922630584 },
{ 0.00000003558, 2.42901552681, 14143.4952424306 },
{ 0.00000003372, 3.86210700128, 1592.5960136328 },
{ 0.00000003374, 0.88776219727, 12036.4607348882 },
{ 0.00000003175, 3.18785710594, 4705.7323075436 },
{ 0.00000003221, 0.61599835472, 8429.2412664666 },
{ 0.00000004132, 5.23992859705, 7084.8967811152 },
{ 0.0000000297, 6.07026318493, 4292.3308329504 },
{ 0.000000029, 2.32464208411, 20.3553193988 },
{ 0.00000003504, 4.79975694359, 6279.5527316424 },
{ 0.0000000295, 1.43108874817, 5746.271337896 },
{ 0.00000002697, 4.80368225199, 7234.794256242 },
{ 0.00000002531, 6.22290682655, 6836.6452528338 },
{ 0.00000002745, 0.93466065396, 5760.4984318976 },
{ 0.0000000325, 3.39954640038, 7632.9432596502 },
{ 0.00000002277, 5.00277837672, 17789.845619785 },
{ 0.00000002075, 3.95534978634, 10213.285546211 },
{ 0.00000002061, 2.22411683077, 5856.4776591154 },
{ 0.00000002252, 5.67166499885, 11499.6562227928 },
{ 0.00000002148, 5.20184578235, 11513.8833167944 },
{ 0.00000001886, 0.53198320577, 3340.6124266998 },
{ 0.00000001875, 4.73511970207, 83996.84731811189 },
{ 0.0000000206, 2.54987293999, 25132.3033999656 },
{ 0.00000001794, 1.47435409831, 4164.311989613 },
{ 0.00000001778, 3.02473091781, 5.5229243074 },
{ 0.00000002029, 0.90960209983, 6256.7775301916 },
{ 0.00000002075, 2.26767270157, 522.5774180938 },
{ 0.00000001772, 3.02622802353, 5753.3848848968 },
{ 0.00000001569, 6.12410242782, 5216.5803728014 },
{ 0.0000000159, 4.63713748247, 3.2863574178 },
{ 0.00000001542, 4.20004448567, 13367.9726311066 },
{ 0.00000001427, 1.19088061711, 3894.1818295422 },
{ 0.00000001375, 3.09301252193, 135.0650800354 },
{ 0.00000001359, 4.24532506641, 426.598190876 },
{ 0.0000000134, 5.76511818622, 6040.3472460174 },
{ 0.00000001284, 3.08524663344, 5643.1785636774 },
{ 0.0000000125, 3.07748157144, 11926.2544136688 },
{ 0.00000001551, 3.07665451458, 6681.2248533996 },
{ 0.00000001268, 2.09196018331, 6290.1893969922 },
{ 0.00000001144, 3.24444699514, 12168.0026965746 },
{ 0.00000001248, 3.44504937285, 536.8045120954 },
{ 0.00000001118, 2.31829670425, 16730.4636895958 },
{ 0.00000001105, 5.31966001019, 23.8784377478 },
{ 0.00000001051, 3.75015946014, 7860.4193924392 },
{ 0.00000001025, 2.44688534235, 1990.745017041 },
{ 0.00000000962, 0.81771017882, 3.881335358 },
{ 0.0000000091, 0.41727865299, 7079.3738568078 },
{ 0.00000000883, 5.16833917651, 11790.6290886588 },
{ 0.00000000957, 4.07673573735, 6127.6554505572 },
{ 0.0000000111, 3.90096793825, 11506.7697697936 },
{ 0.00000000802, 3.88778875582, 10973.55568635 },
{ 0.0000000078, 2.39934293755, 1589.0728952838 },
{ 0.00000000758, 1.30034364248, 103.0927742186 },
{ 0.00000000749, 4.962758033, 6496.3749454294 },
{ 0.00000000765, 3.36312388424, 36.0278666774 },
{ 0.00000000915, 5.41543742089, 206.1855484372 },
{ 0.00000000776, 2.57589093871, 11371.7046897582 },
{ 0.00000000772, 3.98369209464, 955.5997416086 },
{ 0.00000000749, 5.17890001805, 10969.9652576982 },
{ 0.00000000806, 0.34218864254, 9917.6968745098 },
{ 0.00000000728, 5.20962563787, 38.0276726358 },
{ 0.00000000685, 2.77592961854, 20.7753954924 },
{ 0.00000000636, 4.28242193632, 28.4491874678 },
{ 0.00000000608, 5.63278508906, 10984.1923516998 },
{ 0.00000000704, 5.60738823665, 3738.761430108 },
{ 0.00000000685, 0.38876148682, 15.252471185 },
{ 0.00000000601, 0.73489602442, 419.4846438752 },
{ 0.00000000716, 2.65279791438, 6309.3741697912 },
{ 0.00000000584, 5.54502568227, 17298.1823273262 },
{ 0.0000000065, 1.13379656406, 7058.5984613154 },
{ 0.00000000688, 2.59683891779, 3496.032826134 },
{ 0.00000000485, 0.44467180946, 12352.8526045448 },
{ 0.00000000528, 2.74936967681, 3930.2096962196 },
{ 0.00000000597, 5.27668281777, 10575.4066829418 },
{ 0.00000000583, 3.1892906781, 4732.0306273434 },
{ 0.00000000526, 5.01697321546, 5884.9268465832 },
{ 0.0000000054, 1.29175137075, 640.8776073822 },
{ 0.00000000473, 5.4995330697, 5230.807466803 },
{ 0.00000000406, 5.21248452189, 220.4126424388 },
{ 0.00000000395, 1.87474483222, 16200.7727245012 },
{ 0.0000000037, 3.84921354713, 18073.7049386502 },
{ 0.00000000367, 0.88533542778, 6283.14316029419 },
{ 0.00000000379, 0.37983009325, 10177.2576795336 },
{ 0.00000000356, 3.84145204913, 11712.9553182308 },
{ 0.00000000374, 5.01577520608, 7.046236698 },
{ 0.00000000381, 4.30250406634, 6062.6632075526 },
{ 0.00000000471, 0.86381834647, 6069.7767545534 },
{ 0.00000000367, 1.32943839763, 6283.0085396886 },
{ 0.0000000046, 5.19667219575, 6284.0561710596 },
{ 0.00000000333, 5.54256205741, 4686.8894077068 },
{ 0.00000000341, 4.36522989934, 7238.6755916 },
{ 0.00000000336, 4.00205876835, 3097.88382272579 },
{ 0.00000000359, 6.22679790284, 245.8316462294 },
{ 0.00000000307, 2.35299010924, 170.6728706192 },
{ 0.00000000343, 3.77164927143, 6076.8903015542 },
{ 0.00000000296, 5.44152227481, 17260.1546546904 },
{ 0.00000000328, 0.13837875384, 11015.1064773348 },
{ 0.00000000268, 1.1390455063, 12569.6748183318 },
{ 0.00000000263, 0.00538633678, 4136.9104335162 },
{ 0.00000000282, 5.0439983748, 7477.522860216 },
{ 0.00000000288, 3.13401177517, 12559.038152982 },
{ 0.00000000259, 0.93882269387, 5642.1982426092 },
{ 0.00000000292, 1.98420020514, 12132.439962106 },
{ 0.00000000247, 3.84244798532, 5429.8794682394 },
{ 0.00000000245, 5.70467521726, 65147.6197681377 },
{ 0.00000000241, 0.99480969552, 3634.6210245184 },
{ 0.00000000246, 3.06168069935, 110.2063212194 },
{ 0.00000000239, 6.11855909114, 11856.2186514245 },
{ 0.00000000263, 0.66348415419, 21228.3920235458 },
{ 0.00000000262, 1.51070507866, 12146.6670561076 },
{ 0.0000000023, 1.75927314884, 9779.1086761254 },
{ 0.00000000223, 2.00967043606, 6172.869528772 },
{ 0.00000000246, 1.10411690865, 6282.0955289232 },
{ 0.00000000221, 3.03945240854, 8635.9420037632 },
{ 0.00000000214, 4.03840869663, 14314.1681130498 },
{ 0.00000000236, 5.4691507058, 13916.0191096416 },
{ 0.00000000224, 4.68408089456, 24072.9214697764 },
{ 0.00000000212, 2.13695625494, 5849.3641121146 },
{ 0.00000000207, 3.07724246401, 11.729352836 },
{ 0.00000000207, 6.10306282747, 23543.23050468179 },
{ 0.00000000266, 1.00709566823, 2388.8940204492 },
{ 0.00000000217, 6.27837036335, 17267.26820169119 },
{ 0.00000000204, 2.34615348695, 266.6070417218 },
{ 0.00000000195, 5.55015549753, 6133.5126528568 },
};
static double earth_L2[93][3] = {
{ 0.0005291887, 0, 0 },
{ 0.00008719837, 1.07209665242, 6283.0758499914 },
{ 0.00000309125, 0.86728818832, 12566.1516999828 },
{ 0.00000027339, 0.05297871691, 3.523118349 },
{ 0.00000016334, 5.18826691036, 26.2983197998 },
{ 0.00000015752, 3.6845788943, 155.4203994342 },
{ 0.00000009541, 0.75742297675, 18849.2275499742 },
{ 0.00000008937, 2.05705419118, 77713.7714681205 },
{ 0.00000006952, 0.8267330541, 775.522611324 },
{ 0.00000005064, 4.66284525271, 1577.3435424478 },
{ 0.00000004061, 1.03057162962, 7.1135470008 },
{ 0.00000003463, 5.14074632811, 796.2980068164 },
{ 0.00000003169, 6.05291851171, 5507.5532386674 },
{ 0.0000000302, 1.19246506441, 242.728603974 },
{ 0.00000002886, 6.11652627155, 529.6909650946 },
{ 0.0000000381, 3.4405080349, 5573.1428014331 },
{ 0.00000002714, 0.30637881025, 398.1490034082 },
{ 0.00000002371, 4.38118838167, 5223.6939198022 },
{ 0.00000002538, 2.27992810679, 553.5694028424 },
{ 0.00000002079, 3.75435330484, 0.9803210682 },
{ 0.00000001675, 0.90216407959, 951.7184062506 },
{ 0.00000001534, 5.75900462759, 1349.8674096588 },
{ 0.00000001224, 2.97328088405, 2146.1654164752 },
{ 0.00000001449, 4.3641591397, 1748.016413067 },
{ 0.00000001341, 3.72061130861, 1194.4470102246 },
{ 0.00000001254, 2.94846826628, 6438.4962494256 },
{ 0.00000000999, 5.98640014468, 6286.5989683404 },
{ 0.00000000917, 4.79788687522, 5088.6288397668 },
{ 0.00000000828, 3.31321076572, 213.299095438 },
{ 0.00000001103, 1.27104454479, 161000.6857376741 },
{ 0.00000000762, 3.41582762988, 5486.777843175 },
{ 0.00000001044, 0.60409577691, 3154.6870848956 },
{ 0.00000000887, 5.23465144638, 7084.8967811152 },
{ 0.00000000645, 1.60096192515, 2544.3144198834 },
{ 0.00000000681, 3.43155669169, 4694.0029547076 },
{ 0.00000000605, 2.47806340546, 10977.078804699 },
{ 0.00000000706, 6.19393222575, 4690.4798363586 },
{ 0.00000000643, 1.98042503148, 801.8209311238 },
{ 0.00000000502, 1.44394375363, 6836.6452528338 },
{ 0.0000000049, 2.34129524194, 1592.5960136328 },
{ 0.00000000458, 1.30876448575, 4292.3308329504 },
{ 0.00000000431, 0.03526421494, 7234.794256242 },
{ 0.00000000379, 3.17030522615, 6309.3741697912 },
{ 0.00000000348, 0.99049550009, 6040.3472460174 },
{ 0.00000000386, 1.57019797263, 71430.69561812909 },
{ 0.00000000347, 0.67013291338, 1059.3819301892 },
{ 0.00000000458, 3.81499443681, 149854.40013480789 },
{ 0.00000000302, 1.91760044838, 10447.3878396044 },
{ 0.00000000307, 3.55343347416, 8031.0922630584 },
{ 0.00000000395, 4.93701776616, 7632.9432596502 },
{ 0.00000000314, 3.18093696547, 2352.8661537718 },
{ 0.00000000282, 4.41936437052, 9437.762934887 },
{ 0.00000000276, 2.71314254553, 3894.1818295422 },
{ 0.00000000298, 2.5203747421, 6127.6554505572 },
{ 0.0000000023, 1.37790215549, 4705.7323075436 },
{ 0.00000000252, 0.55330133471, 6279.5527316424 },
{ 0.00000000255, 5.26570187369, 6812.766815086 },
{ 0.00000000275, 0.67264264272, 25132.3033999656 },
{ 0.00000000178, 0.92820785174, 1990.745017041 },
{ 0.00000000221, 0.63897368842, 6256.7775301916 },
{ 0.00000000155, 0.77319790838, 14143.4952424306 },
{ 0.0000000015, 2.40470465561, 426.598190876 },
{ 0.00000000196, 6.06877865012, 640.8776073822 },
{ 0.00000000137, 2.21679460145, 8429.2412664666 },
{ 0.00000000127, 3.26094223174, 17789.845619785 },
{ 0.00000000128, 5.47237279946, 12036.4607348882 },
{ 0.00000000122, 2.16291082757, 10213.285546211 },
{ 0.00000000118, 0.45789822268, 7058.5984613154 },
{ 0.00000000141, 2.34932647403, 11506.7697697936 },
{ 0.000000001, 0.85621569847, 6290.1893969922 },
{ 0.00000000092, 5.10587476002, 7079.3738568078 },
{ 0.00000000126, 2.65428307012, 88860.05707098669 },
{ 0.00000000106, 5.85646710022, 7860.4193924392 },
{ 0.00000000084, 3.57457554262, 16730.4636895958 },
{ 0.00000000089, 4.21433259618, 83996.84731811189 },
{ 0.00000000097, 5.57938280855, 13367.9726311066 },
{ 0.00000000102, 2.05853060226, 87.30820453981 },
{ 0.0000000008, 4.73792651816, 11926.2544136688 },
{ 0.0000000008, 5.41418965044, 10973.55568635 },
{ 0.00000000106, 4.10978997399, 3496.032826134 },
{ 0.00000000102, 3.62650006043, 244287.60000722769 },
{ 0.00000000075, 4.89483161769, 5643.1785636774 },
{ 0.00000000087, 0.42863750683, 11015.1064773348 },
{ 0.00000000069, 1.8890876072, 10177.2576795336 },
{ 0.00000000089, 1.35567273119, 6681.2248533996 },
{ 0.00000000066, 0.99455837265, 6525.8044539654 },
{ 0.00000000067, 5.5124099707, 3097.88382272579 },
{ 0.00000000076, 2.72016814799, 4164.311989613 },
{ 0.00000000063, 1.4434990254, 9917.6968745098 },
{ 0.00000000078, 3.51469733747, 11856.2186514245 },
{ 0.00000000085, 0.50956043858, 10575.4066829418 },
{ 0.00000000067, 3.62043033405, 16496.3613962024 },
{ 0.00000000055, 5.24637517308, 3340.6124266998 },
};
static double earth_L3[8][3] = {
{ 0.00000289226, 5.84384198723, 6283.0758499914 },
{ 0.00000034955, 0, 0 },
{ 0.00000016819, 5.48766912348, 12566.1516999828 },
{ 0.00000002962, 5.19577265202, 155.4203994342 },
{ 0.00000001288, 4.72200252235, 3.523118349 },
{ 0.00000000635, 5.96925937141, 242.728603974 },
{ 0.00000000714, 5.30045809128, 18849.2275499742 },
{ 0.00000000402, 3.78682982419, 553.5694028424 },
};
static double earth_L4[4][3] = {
{ 0.00000114084, 3.14159265359, 0 },
{ 0.00000007717, 4.13446589358, 6283.0758499914 },
{ 0.00000000765, 3.83803776214, 12566.1516999828 },
{ 0.0000000042, 0.41925861858, 155.4203994342 },
};
static double earth_L5[4][3] = {
{ 0.00000000878, 3.14159265359, 0 },
{ 0.00000000172, 2.7657906951, 6283.0758499914 },
{ 0.0000000005, 2.01353298182, 155.4203994342 },
{ 0.00000000028, 2.21496423926, 12566.1516999828 },
};
/*
static double earth_B0[20][3] = {
{ 0.0000027962, 3.19870156017, 84334.66158130829 },
{ 0.00000101643, 5.42248619256, 5507.5532386674 },
{ 0.00000080445, 3.88013204458, 5223.6939198022 },
{ 0.00000043806, 3.70444689758, 2352.8661537718 },
{ 0.00000031933, 4.00026369781, 1577.3435424478 },
{ 0.00000022724, 3.9847383156, 1047.7473117547 },
{ 0.00000016392, 3.56456119782, 5856.4776591154 },
{ 0.00000018141, 4.98367470263, 6283.0758499914 },
{ 0.00000014443, 3.70275614914, 9437.762934887 },
{ 0.00000014304, 3.41117857525, 10213.285546211 },
{ 0.00000011246, 4.8282069053, 14143.4952424306 },
{ 0.000000109, 2.08574562327, 6812.766815086 },
{ 0.00000009714, 3.47303947752, 4694.0029547076 },
{ 0.00000010367, 4.05663927946, 71092.88135493269 },
{ 0.00000008775, 4.44016515669, 5753.3848848968 },
{ 0.00000008366, 4.9925151218, 7084.8967811152 },
{ 0.00000006921, 4.32559054073, 6275.9623029906 },
{ 0.00000009145, 1.14182646613, 6620.8901131878 },
{ 0.00000007194, 3.60193205752, 529.6909650946 },
{ 0.00000007698, 5.55425745881, 167621.57585086189 },
};
static double earth_B1[12][3] = {
{ 0.0000000903, 3.8972906189, 5507.5532386674 },
{ 0.00000006177, 1.73038850355, 5223.6939198022 },
{ 0.000000038, 5.24404145734, 2352.8661537718 },
{ 0.00000002834, 2.4734503745, 1577.3435424478 },
{ 0.00000001817, 0.41874743765, 6283.0758499914 },
{ 0.00000001499, 1.83320979291, 5856.4776591154 },
{ 0.00000001466, 5.69401926017, 5753.3848848968 },
{ 0.00000001301, 2.18890066314, 9437.762934887 },
{ 0.00000001233, 4.95222451476, 10213.285546211 },
{ 0.00000001021, 0.12866660208, 7860.4193924392 },
{ 0.00000000982, 0.09005453285, 14143.4952424306 },
{ 0.00000000865, 1.73949953555, 3930.2096962196 },
};
static double earth_B2[4][3] = {
{ 0.00000001662, 1.62703209173, 84334.66158130829 },
{ 0.00000000492, 2.41382223971, 1047.7473117547 },
{ 0.00000000344, 2.24353004539, 5507.5532386674 },
{ 0.00000000258, 6.00906896311, 5223.6939198022 },
};
static double earth_B3[1][3] = {
{ 0, 0, 0 },
};
static double earth_B4[1][3] = {
{ 0, 0, 0 },
};
static double earth_R0[72][3] = {
{ 1.00013988799, 0, 0 },
{ 0.01670699626, 3.09846350771, 6283.0758499914 },
{ 0.00013956023, 3.0552460962, 12566.1516999828 },
{ 0.0000308372, 5.19846674381, 77713.7714681205 },
{ 0.00001628461, 1.17387749012, 5753.3848848968 },
{ 0.00001575568, 2.84685245825, 7860.4193924392 },
{ 0.00000924799, 5.45292234084, 11506.7697697936 },
{ 0.00000542444, 4.56409149777, 3930.2096962196 },
{ 0.0000047211, 3.66100022149, 5884.9268465832 },
{ 0.0000032878, 5.89983646482, 5223.6939198022 },
{ 0.00000345983, 0.96368617687, 5507.5532386674 },
{ 0.00000306784, 0.29867139512, 5573.1428014331 },
{ 0.00000174844, 3.01193636534, 18849.2275499742 },
{ 0.00000243189, 4.27349536153, 11790.6290886588 },
{ 0.00000211829, 5.84714540314, 1577.3435424478 },
{ 0.00000185752, 5.02194447178, 10977.078804699 },
{ 0.00000109835, 5.05510636285, 5486.777843175 },
{ 0.00000098316, 0.88681311277, 6069.7767545534 },
{ 0.00000086499, 5.68959778254, 15720.8387848784 },
{ 0.00000085825, 1.27083733351, 161000.6857376741 },
{ 0.00000062916, 0.92177108832, 529.6909650946 },
{ 0.00000057056, 2.01374292014, 83996.84731811189 },
{ 0.00000064903, 0.27250613787, 17260.1546546904 },
{ 0.00000049384, 3.24501240359, 2544.3144198834 },
{ 0.00000055736, 5.24159798933, 71430.69561812909 },
{ 0.00000042515, 6.01110242003, 6275.9623029906 },
{ 0.00000046963, 2.57805070386, 775.522611324 },
{ 0.00000038968, 5.36071738169, 4694.0029547076 },
{ 0.00000044661, 5.53715807302, 9437.762934887 },
{ 0.0000003566, 1.67468058995, 12036.4607348882 },
{ 0.00000031921, 0.18368229781, 5088.6288397668 },
{ 0.00000031846, 1.77775642085, 398.1490034082 },
{ 0.00000033193, 0.24370300098, 7084.8967811152 },
{ 0.00000038245, 2.39255343974, 8827.3902698748 },
{ 0.00000028464, 1.21344868176, 6286.5989683404 },
{ 0.0000003749, 0.82952922332, 19651.048481098 },
{ 0.00000036957, 4.90107591914, 12139.5535091068 },
{ 0.00000034537, 1.84270693282, 2942.4634232916 },
{ 0.00000026275, 4.58896850401, 10447.3878396044 },
{ 0.00000024596, 3.78660875483, 8429.2412664666 },
{ 0.00000023587, 0.26866117066, 796.2980068164 },
{ 0.00000027793, 1.89934330904, 6279.5527316424 },
{ 0.00000023927, 4.99598548138, 5856.4776591154 },
{ 0.00000020349, 4.65267995431, 2146.1654164752 },
{ 0.00000023287, 2.80783650928, 14143.4952424306 },
{ 0.00000022103, 1.95004702988, 3154.6870848956 },
{ 0.00000019506, 5.38227371393, 2352.8661537718 },
{ 0.00000017958, 0.19871379385, 6812.766815086 },
{ 0.00000017174, 4.43315560735, 10213.285546211 },
{ 0.0000001619, 5.23160507859, 17789.845619785 },
{ 0.00000017314, 6.15200787916, 16730.4636895958 },
{ 0.00000013814, 5.18962074032, 8031.0922630584 },
{ 0.00000018833, 0.67306674027, 149854.40013480789 },
{ 0.00000018331, 2.25348733734, 23581.2581773176 },
{ 0.00000013641, 3.68516118804, 4705.7323075436 },
{ 0.00000013139, 0.65289581324, 13367.9726311066 },
{ 0.00000010414, 4.33285688538, 11769.8536931664 },
{ 0.00000009978, 4.20126336355, 6309.3741697912 },
{ 0.00000010169, 1.59390681369, 4690.4798363586 },
{ 0.00000007564, 2.6256059739, 6256.7775301916 },
{ 0.00000009661, 3.6758679122, 27511.4678735372 },
{ 0.00000006743, 0.56270332741, 3340.6124266998 },
{ 0.00000008743, 6.06359123461, 1748.016413067 },
{ 0.00000007786, 3.67371235637, 12168.0026965746 },
{ 0.00000006633, 5.66149277792, 11371.7046897582 },
{ 0.00000007712, 0.31242577789, 7632.9432596502 },
{ 0.00000006592, 3.13576266188, 801.8209311238 },
{ 0.0000000746, 5.64757188143, 11926.2544136688 },
{ 0.00000006933, 2.923845864, 6681.2248533996 },
{ 0.00000006802, 1.4232980642, 23013.5395395872 },
{ 0.00000006115, 5.13393615454, 1194.4470102246 },
{ 0.00000006477, 2.64986648492, 19804.8272915828 },
};
static double earth_R1[10][3] = {
{ 0.00103018608, 1.10748969588, 6283.0758499914 },
{ 0.00001721238, 1.06442301418, 12566.1516999828 },
{ 0.00000702215, 3.14159265359, 0 },
{ 0.00000032346, 1.02169059149, 18849.2275499742 },
{ 0.00000030799, 2.84353804832, 5507.5532386674 },
{ 0.00000024971, 1.31906709482, 5223.6939198022 },
{ 0.00000018485, 1.42429748614, 1577.3435424478 },
{ 0.00000010078, 5.91378194648, 10977.078804699 },
{ 0.00000008634, 0.27146150602, 5486.777843175 },
{ 0.00000008654, 1.42046854427, 6275.9623029906 },
};
static double earth_R2[6][3] = {
{ 0.00004359385, 5.78455133738, 6283.0758499914 },
{ 0.00000123633, 5.57934722157, 12566.1516999828 },
{ 0.00000012341, 3.14159265359, 0 },
{ 0.00000008792, 3.62777733395, 77713.7714681205 },
{ 0.00000005689, 1.86958905084, 5573.1428014331 },
{ 0.00000003301, 5.47027913302, 18849.2275499742 },
};
static double earth_R3[2][3] = {
{ 0.00000144595, 4.27319435148, 6283.0758499914 },
{ 0.00000006729, 3.91697608662, 12566.1516999828 },
};
static double earth_R4[1][3] = {
{ 0.00000003858, 2.56384387339, 6283.0758499914 },
};
static double earth_R5[1][3] = {
{ 0, 0, 0 },
};
*/
/* helper function for calculate VSOP87 */
double vsopLx(double vsopterms[][3], size_t rowcount, double t)
{
double lx = 0;
int i;
for (i = 0; i < rowcount; i++) {
lx += vsopterms[i][0] * cos(vsopterms[i][1] + vsopterms[i][2] * t);
}
return lx;
}
/* Calculate ecliptical longitude of earth in heliocentric coordinates,
* use VSOP87D table, heliocentric spherical, coordinates referred to the mean
* equinox of the date,
*
* In A&A, Meeus said while the complete VSOP87 yields an accuracy of 0.01",
* the A&A abridge VSOP87 has an accuracy of 1" for -2000 - +6000.
*
* The VSOP87D table used here is a truncated version, done by the
* vsoptrunc-sph.c from Celestia.
*
* Arg:
* jd: in JDTT
* Return:
* earth longitude in radians, referred to mean dynamical ecliptic and
* equinox of the date
*/
double vsop(double jd)
{
double t, lon;
double L0, L1, L2, L3, L4, L5;
t = (jd - J2000) / 365250.0;
L0 = vsopLx(earth_L0, sizeof(earth_L0) / 24, t);
L1 = vsopLx(earth_L1, sizeof(earth_L1) / 24, t);
L2 = vsopLx(earth_L2, sizeof(earth_L2) / 24, t);
L3 = vsopLx(earth_L3, sizeof(earth_L3) / 24, t);
L4 = vsopLx(earth_L4, sizeof(earth_L4) / 24, t);
L5 = vsopLx(earth_L5, sizeof(earth_L5) / 24, t);
lon = (L0 + t * (L1 + t * (L2 + t * (L3 + t * (L4 + t * L5)))));
/* adjust FK5 */
lon += -4.379321981462438e-07;
return lon;
}
/* calculate the apprent place of the Sun.
* Arg:
* jd as jd
* Return:
* geocentric longitude in radians
*/
double apparentsun(double jd, int ignorenutation)
{
double geolon;
geolon = vsop(jd) + PI;
/* compensate nutation */
if (!ignorenutation)
geolon += nutation(jd);
geolon += lightabbr_high(jd);
return geolon;
}