-
Notifications
You must be signed in to change notification settings - Fork 73
/
Copy pathmoon_north_south.js
208 lines (174 loc) · 6.57 KB
/
moon_north_south.js
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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
/*
moon_north_south.js - by Don Cross - 2023-08-22
Calculates when the Moon reaches extreme declination
(most north or south with respect to the plane of the Earth's equator)
and when the Moon reaches the most extreme ecliptic latitude
(most north or south with respect to the plane of the Earth's orbit around the Sun).
*/
const Astronomy = require('./astronomy.js');
function Ecliptic(time) {
// Return the Moon's ecliptic latitude at the given time.
const vec = Astronomy.GeoMoon(time);
const ecl = Astronomy.Ecliptic(vec);
return ecl.elat;
}
function Equatorial(time) {
// Return the Moon's declination angle at the given time.
// Start with the Moon's position vector in J2000 coordinates.
const eqj = Astronomy.GeoMoon(time);
// Find rotation matrix to convert J2000 coordinates to equator-of-date.
const rot = Astronomy.Rotation_EQJ_EQD(time);
// Transform coordinates into equator-of-date.
const eqd = Astronomy.RotateVector(rot, eqj);
// Convert to angular coordinates to find declination angle.
const equ = Astronomy.EquatorFromVector(eqd);
return equ.dec;
}
function Search(startTime, direction, func) {
// Create a callback function that reports the rate of change of the desired variable.
function f(t) {
const dt = 1 / 86400; // one second, expressed in days
const x1 = func(t.AddDays(-dt));
const x2 = func(t.AddDays(+dt));
return direction*(x2 - x1);
}
// Search forward 10 days at a time until we find a solution.
// Because the Moon's orbit takes about 29 days, we want an interval
// that is less than half that amount of time. This prevents
// finding more than one extreme (minimum/maximum) in a single
// search interval, which would cause the search to fail.
let t1 = startTime;
while (true) {
const t2 = t1.AddDays(10);
const tx = Astronomy.Search(f, t1, t2);
if (tx) {
return tx; // found a solution!
}
t1 = t2;
}
}
function Solve(time1, direction, func, comment) {
const time = Search(time1, direction, func);
const angle = func(time);
console.log(`${time} Moon next reaches ${comment} = ${angle.toFixed(7)}.`);
}
function ParseDate(text) {
const d = new Date(text);
if (!Number.isFinite(d.getTime())) {
console.error(`ERROR: Not a valid date: "${text}"`);
process.exit(1);
}
return d;
}
function Demo() {
if (process.argv.length < 3) {
console.log('USAGE: node moon_north_south yyyy-mm-dd[Thh:mm:ssZ] ...');
process.exit(1);
} else {
for (let i = 2; i < process.argv.length; ++i) {
const time1 = Astronomy.MakeTime(ParseDate(process.argv[i]));
console.log(`${time1} Starting search.`);
Solve(time1, -1, Ecliptic, 'maximum ecliptic latitude');
Solve(time1, +1, Ecliptic, 'minimum ecliptic latitude');
Solve(time1, -1, Equatorial, 'maximum declination');
Solve(time1, +1, Equatorial, 'minimum declination');
console.log();
}
process.exit(0);
}
}
Demo();
/*
---------------------------------------------------------------------------------------
JPL Horizons data for maximum ecliptic latitude:
time longitude latitude
2023-Sep-10 07:35 115.5157583 5.1662606
2023-Sep-10 07:40 115.5572561 5.1662659
2023-Sep-10 07:45 115.5987527 5.1662686 <== max
2023-Sep-10 07:50 115.6402480 5.1662685
2023-Sep-10 07:55 115.6817421 5.1662657
This program says:
2023-09-10T07:47:00.420Z Moon next reaches maximum ecliptic latitude = 5.1662436.
---------------------------------------------------------------------------------------
JPL Horizons data for minimum ecliptic latitude:
time longitude latitude
2023-Aug-28 08:15 296.0963021 -5.1108070
2023-Aug-28 08:20 296.1479665 -5.1108143
2023-Aug-28 08:25 296.1996351 -5.1108175 <== min
2023-Aug-28 08:30 296.2513079 -5.1108166
2023-Aug-28 08:35 296.3029851 -5.1108115
This program says:
2023-08-28T08:26:14.233Z Moon next reaches minimum ecliptic latitude = -5.1107593.
---------------------------------------------------------------------------------------
JPL Horizons data for maximum declination:
time RA(deg) DEC(deg)
2023-Sep-08 13:00 94.62655 28.17823
2023-Sep-08 13:05 94.67459 28.17827
2023-Sep-08 13:10 94.72262 28.17829
2023-Sep-08 13:15 94.77064 28.17829 <== max
2023-Sep-08 13:20 94.81867 28.17828
2023-Sep-08 13:25 94.86669 28.17824
This program says:
2023-09-08T13:13:01.641Z Moon next reaches maximum declination = 28.1783302.
---------------------------------------------------------------------------------------
JPL Horizons data for maximum declination:
time RA(deg) DEC(deg)
2023-Aug-26 20:05 274.66967 -28.10661
2023-Aug-26 20:10 274.72558 -28.10665
2023-Aug-26 20:15 274.78150 -28.10668
2023-Aug-26 20:20 274.83743 -28.10668 <== min
2023-Aug-26 20:25 274.89336 -28.10667
2023-Aug-26 20:30 274.94930 -28.10662
This program says:
2023-08-26T20:18:11.883Z Moon next reaches minimum declination = -28.1066748.
---------------------------------------------------------------------------------------
JPL Batch Data for calculating ecliptic latitudes:
!$$SOF
MAKE_EPHEM=YES
COMMAND=301
EPHEM_TYPE=OBSERVER
CENTER='500@399'
START_TIME='2023-08-28'
STOP_TIME='2023-08-29'
STEP_SIZE='5 MINUTES'
QUANTITIES='31'
REF_SYSTEM='ICRF'
CAL_FORMAT='CAL'
CAL_TYPE='M'
TIME_DIGITS='MINUTES'
ANG_FORMAT='HMS'
APPARENT='AIRLESS'
RANGE_UNITS='AU'
SUPPRESS_RANGE_RATE='NO'
SKIP_DAYLT='NO'
SOLAR_ELONG='0,180'
EXTRA_PREC='NO'
R_T_S_ONLY='NO'
CSV_FORMAT='NO'
OBJ_DATA='YES'
---------------------------------------------------------------------------------------
JPL Batch data for calculating equatorial coordinates:
!$$SOF
MAKE_EPHEM=YES
COMMAND=301
EPHEM_TYPE=OBSERVER
CENTER='500@399'
START_TIME='2023-08-26'
STOP_TIME='2023-08-27'
STEP_SIZE='5 MINUTES'
QUANTITIES='2'
REF_SYSTEM='ICRF'
CAL_FORMAT='CAL'
CAL_TYPE='M'
TIME_DIGITS='MINUTES'
ANG_FORMAT='DEG'
APPARENT='AIRLESS'
RANGE_UNITS='AU'
SUPPRESS_RANGE_RATE='NO'
SKIP_DAYLT='NO'
SOLAR_ELONG='0,180'
EXTRA_PREC='NO'
R_T_S_ONLY='NO'
CSV_FORMAT='NO'
OBJ_DATA='YES'
*/