1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
|
#include "project.h"
#define D2R(a) ((a)*(M_PI/180.))
// shamelessly stolen from Meeus Astronmical Algorithms Chapter 27 table 27.B
// chapter 25 is probably more tractable, but this is easier to code up
# if 0
static double mean_vernal_equinox (unsigned year)
{
double y = year;
y -= 2000;
y *= 0.001;
return 2451623.80984 + 365242.37404 * y + 0.05169 * (y * y) - 0.00411 * (y * y * y) - 0.00057 * (y * y * y * y);
}
static double mean_summer_solstice (unsigned year)
{
double y = year;
y -= 2000;
y *= 0.001;
return 2451716.56767 + 365241.62603 * y + 0.00325 * (y * y) + 0.00888 * (y * y * y) - 0.00030 * (y * y * y * y);
}
#endif
static double mean_autumnal_equinox (unsigned year)
{
double y = year;
y -= 2000;
y *= 0.001;
return 2451810.21715 + 365242.01767 * y - 0.11575 * (y * y) + 0.00337 * (y * y * y) + 0.00078 * (y * y * y * y);
}
#if 0
static double mean_winter_solstice (unsigned year)
{
double y = year;
y -= 2000;
y *= 0.001;
return 2451900.05952 + 365242.74049 * y - 0.06223 * (y * y) - 0.00823 * (y * y * y) + 0.00032 * (y * y * y * y);
}
#endif
static double orbital_periodic_terms (double t)
{
#define N 24
const double A[N] = {485, 203, 199, 182, 156, 136, 77, 74, 70, 58, 52, 50, 45,
44, 29, 18, 17, 16, 14, 12, 12, 12, 9, 8
};
const double B[N] = {D2R (324.96), D2R (337.23), D2R (342.08), D2R (27.85),
D2R (73.14), D2R (171.52), D2R (222.54), D2R (296.72),
D2R (243.58), D2R (119.81), D2R (297.17), D2R (21.02),
D2R (247.54), D2R (325.15), D2R (60.93), D2R (155.12),
D2R (288.79), D2R (198.04), D2R (199.76), D2R (95.39),
D2R (287.11), D2R (320.81), D2R (227.73), D2R (15.45)
};
const double C[N] = {D2R (1934.136), D2R (32964.467), D2R (20.186),
D2R (445267.112), D2R (45036.886), D2R (22518.443),
D2R (65928.934), D2R (3034.906), D2R (9037.513),
D2R (33718.147), D2R (150.678), D2R (2281.226),
D2R (29929.562), D2R (31555.956), D2R (4443.417),
D2R (67555.328), D2R (4562.452), D2R (62894.029),
D2R (31436.921), D2R (14577.848), D2R (31931.756),
D2R (34777.259), D2R (1222.114), D2R (16859.074)
};
double s = 0;
unsigned i;
for (i = 0; i < N; ++i)
s += A[i] * cos (B[i] + (C[i] * t));
return s;
}
static double mean_to_real (double j0)
{
double t = (j0 - 2451545.) / 36525.;
double w = D2R ((35999.373 * t) - 2.47);
double dl = 1 + 0.0334 * cos (w) + 0.0007 * cos (2. * w);
double s = orbital_periodic_terms (t);
#if 0
printf ("j0=%.6f\r\n", j0);
printf ("t=%.6f\r\n", t);
printf ("w=%.6f\r\n", w);
printf ("dl=%.6f\r\n", dl);
printf ("s=%.6f\r\n", s);
#endif
return j0 + ((0.00001 * s) / dl);
}
double autumnal_equinox (unsigned y)
{
return mean_to_real (mean_autumnal_equinox (y));
// return mean_to_real (mean_summer_solstice (y));
}
|
