-
Notifications
You must be signed in to change notification settings - Fork 0
/
planets_with_keplers_equation.html
376 lines (321 loc) · 13.9 KB
/
planets_with_keplers_equation.html
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
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
<html>
<head>
<title>Compute Planet Positions Using Kepler's Equation</title>
<link rel="stylesheet" href="default.css">
<link rel="stylesheet" href="highlight/styles/default.css">
<meta name="viewport" content="width=device-width, initial-scale=1" />
</head>
<body>
<h1>Compute Planet Positions Using Kepler's Equation</h1>
<p>This is an implementation of the algorithm described in Chapter 8 of the Explanatory Supplement to the Astronomical Almanac 3rd ed P340.
Two sets of elements are provided, and they are accurate for different intervals. The first set is valid for 1800-2050, and the worst
case accuracy is for Saturn at about 10 arc minutes. Most others are under 1 arcmin. The second set is valid for 3000BC to 3000AD, and has
a worst case accuracy (for Uranus) of about half a degree. Most others are on the order of 1 arc minute. See table 8.10 in the Explanatory
Supplement for full accuracy details.</p>
<p>You'll notice that each planet has 12 elements rather than the nomral 6 Kelperian elements. A pure Kelperian element solve looses accuracy
rather quickly due to purturation effects from other planets. So the authors of this routine have provided an additional 6 elements which are
corrections to the Kelperian elements based on time since Jan 1 2000. The corrections were specifically fit to the specified time periods,
so they are considered invalid outside of that time period.
</p>
<p>An explanation of the method used is also available as
<a href="https://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf">Keperian Elements for Approximat Positions of the Major Planets</a>.
A more general overview is at <a href="https://ssd.jpl.nasa.gov/?planet_pos">SSD JPL</a>.</p>
<p>Since this routine contains only elements for the Earth-Moon Barycenter, and not the Earth or the Moon, <b>this example
uses the Earth-Moon Barycenter as the location for the observer</b>. Due to parallax, the positions will be off slightly from those
generated for a Geocentric or Topocentric observer. Since the planets are far away, this error will be insignificant for many purposes.
</p>
<p>Due to it's compactness, and relatively quick computation, this method can be quite useful in a lot of applications where the high
memory footprint, and number of computations of other ephemeris methods are impractical or unnecessary.
</p>
<form>
<input type="datetime-local" onchange="onDateChanged()" id="inputDate" step="1"> Local time
</form>
<div id="outputJulianDate"></div>
<table border="1" cellspacing="0" cellpadding="5" id="outputTable">
</table>
<pre id="sourceOutput">
<code class="JavaScript" id='code1'></code>
</pre>
<script id='sourceCode'>
/*
Greg Miller gmiller@gregmiller.net 2021
http://www.celestialprogramming.com/
Released as public domain
*/
const data1800to2050 = [ //https://ssd.jpl.nasa.gov/txt/p_elem_t1.txt
//Mercury
[[ 0.38709927, 0.20563593, 7.00497902, 252.25032350, 77.45779628, 48.33076593],
[ 0.00000037, 0.00001906, -0.00594749, 149472.67411175, 0.16047689, -0.12534081]],
//Venus,
[[ 0.72333566, 0.00677672, 3.39467605, 181.97909950, 131.60246718, 76.67984255],
[ 0.00000390, -0.00004107, -0.00078890, 58517.81538729, 0.00268329, -0.27769418]],
//EM Bary,
[[ 1.00000261, 0.01671123, -0.00001531, 100.46457166, 102.93768193, 0.0],
[ 0.00000562, -0.00004392, -0.01294668, 35999.37244981, 0.32327364, 0.0]],
//Mars,
[[ 1.52371034, 0.09339410, 1.84969142, -4.55343205, -23.94362959, 49.55953891],
[ 0.00001847, 0.00007882, -0.00813131, 19140.30268499, 0.44441088, -0.29257343]],
//Jupiter,
[[ 5.20288700, 0.04838624, 1.30439695, 34.39644051, 14.72847983, 100.47390909],
[-0.00011607, -0.00013253, -0.00183714, 3034.74612775, 0.21252668, 0.20469106]],
//Saturn,
[[ 9.53667594, 0.05386179, 2.48599187, 49.95424423, 92.59887831, 113.66242448],
[-0.00125060, -0.00050991, 0.00193609, 1222.49362201, -0.41897216, -0.28867794]],
//Uranus,
[[19.18916464, 0.04725744, 0.77263783, 313.23810451, 170.95427630, 74.01692503],
[-0.00196176, -0.00004397, -0.00242939, 428.48202785, 0.40805281, 0.04240589]],
//Neptune,
[[30.06992276, 0.00859048, 1.77004347, -55.12002969, 44.96476227, 131.78422574],
[ 0.00026291, 0.00005105, 0.00035372, 218.45945325, -0.32241464, -0.00508664]],
//Pluto,
[[39.48211675, 0.24882730, 17.14001206, 238.92903833, 224.06891629, 110.30393684],
[-0.00031596, 0.00005170, 0.00004818, 145.20780515, -0.04062942, -0.01183482]]
];
const data3000BCto3000AD = [ //https://ssd.jpl.nasa.gov/txt/p_elem_t2.txt
//Mercury
[[ 0.38709843, 0.20563661, 7.00559432, 252.25166724, 77.45771895, 48.33961819],
[ 0.00000000, 0.00002123, -0.00590158, 149472.67486623, 0.15940013, -0.12214182]],
//Venus
[[ 0.72332102, 0.00676399, 3.39777545, 181.97970850, 131.76755713, 76.67261496],
[-0.00000026, -0.00005107, 0.00043494, 58517.81560260, 0.05679648, -0.27274174]],
//EM Bary
[[ 1.00000018, 0.01673163, -0.00054346, 100.46691572, 102.93005885, -5.11260389],
[-0.00000003, -0.00003661, -0.01337178, 35999.37306329, 0.31795260, -0.24123856]],
//Mars
[[ 1.52371243, 0.09336511, 1.85181869, -4.56813164, -23.91744784, 49.71320984],
[ 0.00000097, 0.00009149, -0.00724757, 19140.29934243, 0.45223625, -0.26852431]],
//Jupiter
[[ 5.20248019, 0.04853590, 1.29861416, 34.33479152, 14.27495244, 100.29282654],
[-0.00002864, 0.00018026, -0.00322699, 3034.90371757, 0.18199196, 0.13024619]],
//Saturn
[[ 9.54149883, 0.05550825, 2.49424102, 50.07571329, 92.86136063, 113.63998702],
[-0.00003065, -0.00032044, 0.00451969, 1222.11494724, 0.54179478, -0.25015002]],
//Uranus
[[19.18797948, 0.04685740, 0.77298127, 314.20276625, 172.43404441, 73.96250215],
[-0.00020455, -0.00001550, -0.00180155, 428.49512595, 0.09266985, 0.05739699]],
//Neptune
[[30.06952752, 0.00895439, 1.77005520, 304.22289287, 46.68158724, 131.78635853],
[ 0.00006447, 0.00000818, 0.00022400, 218.46515314, 0.01009938, -0.00606302]],
//Pluto
[[39.48686035, 0.24885238, 17.14104260, 238.96535011, 224.09702598, 110.30167986],
[ 0.00449751, 0.00006016, 0.00000501, 145.18042903, -0.00968827, -0.00809981]],
];
// b c s f
const extraTerms = [
[-0.00012452, 0.06064060, -0.35635438, 38.35125000], //Jupiter
[ 0.00025899, -0.13434469, 0.87320147, 38.35125000], //Saturn
[ 0.00058331, -0.97731848, 0.17689245, 7.67025000], //Uranus
[-0.00041348, 0.68346318, -0.10162547, 7.67025000], //Neptune
[-0.01262724, 0, 0, 0,] //Pluto
];
const planetNames=["Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune","Pluto","Sun"];
function computePlanetPosition(jd,elements,rates,extraTerms){
const toRad=Math.PI/180;
//Algorithm from Explanatory Supplement to the Astronomical Almanac ch8 P340
//Step 1:
const T=(jd-2451545.0)/36525;
let a=elements[0]+rates[0]*T;
let e=elements[1]+rates[1]*T;
let I=elements[2]+rates[2]*T;
let L=elements[3]+rates[3]*T;
let w=elements[4]+rates[4]*T;
let O=elements[5]+rates[5]*T;
//Step 2:
let ww=w-O;
let M=L - w;
if(extraTerms.length > 0){
const b=extraTerms[0];
const c=extraTerms[1];
const s=extraTerms[2];
const f=extraTerms[3];
M=L - w + b*T*T + c*Math.cos(f*T*toRad) + s*Math.sin(f*T*toRad);
}
while(M>180){M-=360;}
let E=M+57.29578*e*Math.sin(M*toRad);
let dE=1;
let n=0;
while(dE>1e-7 && n<10){
dE=SolveKepler(M,e,E);
E+=dE;
n++;
}
//Step 4:
const xp=a*(Math.cos(E*toRad)-e);
const yp=a*Math.sqrt(1-e*e)*Math.sin(E*toRad);
const zp=0;
//Step 5:
a*=toRad; e*=toRad; I*=toRad; L*=toRad; ww*=toRad; O*=toRad;
const xecl=(Math.cos(ww)*Math.cos(O)-Math.sin(ww)*Math.sin(O)*Math.cos(I))*xp + (-Math.sin(ww)*Math.cos(O)-Math.cos(ww)*Math.sin(O)*Math.cos(I))*yp;
const yecl=(Math.cos(ww)*Math.sin(O)+Math.sin(ww)*Math.cos(O)*Math.cos(I))*xp + (-Math.sin(ww)*Math.sin(O)+Math.cos(ww)*Math.cos(O)*Math.cos(I))*yp;
const zecl=(Math.sin(ww)*Math.sin(I))*xp + (Math.cos(ww)*Math.sin(I))*yp;
//Step 6:
const eps=23.43928*toRad;
const x=xecl;
const y=Math.cos(eps)*yecl - Math.sin(eps)*zecl;
const z=Math.sin(eps)*yecl + Math.cos(eps)*zecl;
return [x,y,z];
}
function SolveKepler(M,e,E){
const toRad=Math.PI/180;
dM=M - (E-e/toRad*Math.sin(E*toRad));
dE=dM/(1-e*Math.cos(E*toRad));
return dE;
}
</script>
<script>
computeAll(new Date());
function onDateChanged(){
let d=document.getElementById("inputDate").value;
computeAll(new Date(d));
}
function computeAll(d){
const t=d.getTime();
const jd=JulianDateFromUnixTime(t);
document.getElementById("outputJulianDate").innerHTML=d.toUTCString()+"<br>Julian Date: "+jd;
let planetS=new Array();
let planetL=new Array();
for(let i=0;i<data1800to2050.length;i++){
planetS[i]=computePlanetShort(i,jd);
planetL[i]=computePlanetLong(i,jd);
}
clearResultsTable();
for(let i=0;i<data1800to2050.length;i++){
if(i!=2){
planetS[i]=sub(planetS[i],planetS[2]);
planetL[i]=sub(planetL[i],planetL[2]);
printPlanet(planetS[i],i,false);
printPlanet(planetL[i],i,true);
}
}
let sunS=new Array();
sunS[0]=-planetS[2][0];
sunS[1]=-planetS[2][1];
sunS[2]=-planetS[2][2];
printPlanet(sunS,9,false)
let sunL=new Array();
sunL[0]=-planetL[2][0];
sunL[1]=-planetL[2][1];
sunL[2]=-planetL[2][2];
printPlanet(sunL,9,true)
}
function clearResultsTable(){
const t = document.getElementById("outputTable");
while(t.hasChildNodes())
{
t.removeChild(t.firstChild);
}
const r=t.insertRow(t.rows.length);
const c1=r.insertCell(0);
c1.innerHTML="Planet";
const c11=r.insertCell(1);
c11.innerHTML="Series";
const c2=r.insertCell(2);
c2.innerHTML="RA (J2000)";
const c3=r.insertCell(3);
c3.innerHTML="DEC (J2000)";
const c4=r.insertCell(4);
c4.innerHTML="Distance AU";
}
function sub(a,b){
let temp=new Array();
for(let i=0;i<a.length;i++){
temp[i]=a[i]-b[i];
}
return temp;
}
function printPlanet(a,planet,longTerm){
const b=xyzToRaDec(a);
b[0]*=180/Math.PI;
b[1]*=180/Math.PI;
const t=document.getElementById("outputTable");
const r=t.insertRow(t.rows.length);
const c1=r.insertCell(0);
c1.innerHTML=planetNames[planet];
const c11=r.insertCell(1);
c11.innerHTML="Short";
if(longTerm){
c11.innerHTML="Long";
}
const c2=r.insertCell(2);
c2.innerHTML=degreesToHMS(b[0]);
const c3=r.insertCell(3);
c3.innerHTML=degreesToDMS(b[1]);
const c4=r.insertCell(4);
c4.innerHTML=b[2]+" AU";
}
function degreesToDMS(d){
let sign="+";
if(d<0){sign="-";}
d=Math.abs(d);
let deg=Math.floor(d);
d-=deg;
d*=60;
let min=Math.floor(d);
d-=min;
d*=60;
let sec=d.toFixed(2);
if(deg<10){deg="0"+deg;}
if(min<10){min="0"+min;}
if(sec<10){sec="0"+sec;}
return sign+deg+"° "+min+"' "+sec+"\"";
}
function degreesToHMS(d){
d=d/15;
let sign="+";
if(d<0){sign="-";}
d=Math.abs(d);
let deg=Math.floor(d);
d-=deg;
d*=60;
let min=Math.floor(d);
d-=min;
d*=60;
let sec=d.toFixed(2);
if(deg<10){deg="0"+deg;}
if(min<10){min="0"+min;}
if(sec<10){sec="0"+sec;}
return sign+deg+"h "+min+"m "+sec+"s";
}
function computePlanetShort(planet,jd){
return computePlanetPosition(jd,data1800to2050[planet][0],data1800to2050[planet][1],[]);
}
function computePlanetLong(planet,jd){
let extra=[];
if(planet>3){
extra=extraTerms[planet-4];
}
return computePlanetPosition(jd,data3000BCto3000AD[planet][0],data3000BCto3000AD[planet][1],extra);
}
function xyzToRaDec(target){
const x=target[0];
const y=target[1];
const z=target[2];
//Convert from Cartesian to polar coordinates
const r=Math.sqrt(x*x+y*y+z*z);
let l=Math.atan2(y,x);
let t=Math.acos(z/r);
//Make sure RA is positive, and Dec is in range +/-90
if(l<0){l+=2*Math.PI;}
t=.5*Math.PI-t;
return [l,t,r];
}
function JulianDateFromUnixTime(t){
//Not valid for dates before Oct 15, 1582
return (t / 86400000) + 2440587.5;
}
function UnixTimeFromJulianDate(jd){
//Not valid for dates before Oct 15, 1582
return (jd-2440587.5)*86400000;
}
</script>
<script>
let t=document.getElementById("sourceCode").innerText;
t=t.replaceAll("&","&")
t=t.replaceAll("<","<")
t=t.replaceAll(">",">")
document.getElementById("code1").innerHTML=t;
</script>
<script src="highlight/highlight.pack.js"></script>
<script>hljs.initHighlightingOnLoad();</script>
</body>
</html>