02-12-2018, 09:53 AM
From a recent conversation with a friend's on the Kepler's laws - his son's subject for a High school paper. The goal was to calculate the value of the swept area between two instants, as determined by the azimuth angles (a1, a2) of the segments linking the focus of the ellipse with the planet at those moments.
Initially I thought the formulas would involve the Elliptic functions, as elliptical sectors were involved - but it appears that's not the case when the coordinates are centered at the focal point, instead of at the center of the ellipse. I found that fact interesting, as it only involves trigonometric functions (not even hyperbolic).
Here's the reference I followed to program it, a good article that describes an ingenious approach - avoiding painful integration steps. It may not be the simplest way to get this done, chime in if you know a better one.
http://www.bado-shanai.net/Platonic%20Dr...Sector.htm
And here's the FOCAL listing for a plain HP-41 - no extensions whatsoever. The result is the area swept between the two positions defined by the angles a1 and a2; a2 > a1. The parameters a,b are the semi-axis of the ellipse, a>b.
(*) the symbol "<)" is for the angle character.
Example: calculate the area swept between a1 = pi/4 and a2 = 3.pi/4, if the ellipse parameters are a= 2, b= 3
Solution: A = 1.989554087
Once the area is obtained, and knowing the period of the orbiting (T), it's straightforward to determine the time taken by the planet to travel between the two positions, with the direct application of Kepler's 2nd. law:
t = T . Area / pi.a.b
Cheers,
ÁM
Initially I thought the formulas would involve the Elliptic functions, as elliptical sectors were involved - but it appears that's not the case when the coordinates are centered at the focal point, instead of at the center of the ellipse. I found that fact interesting, as it only involves trigonometric functions (not even hyperbolic).
Here's the reference I followed to program it, a good article that describes an ingenious approach - avoiding painful integration steps. It may not be the simplest way to get this done, chime in if you know a better one.
http://www.bado-shanai.net/Platonic%20Dr...Sector.htm
And here's the FOCAL listing for a plain HP-41 - no extensions whatsoever. The result is the area swept between the two positions defined by the angles a1 and a2; a2 > a1. The parameters a,b are the semi-axis of the ellipse, a>b.
Code:
1 LBL "K2+"
2 RAD
3 "a^b=?"
4 PROMPT
5 X>Y? b > a ?
6 X<>Y yes, swap
7 STO 02 b
8 X^2 b^2
9 X<>Y
10 STO 01 a
11 X^2 a^2
12 - a^2 - b^2
13 ABS
14 SQRT
15 RCL 01 a
16 /
17 STO 00 e
18 LBL 00
19 "<)1^<)2=?"
20 PROMPT
21 X<Y? a2 < a1 ?
22 X<>Y yes, swap
23 STO 04 a2
24 RDN
25 STO 03 a1
26 PI
27 X<=Y? a1 >= pi ?
28 GTO 01 yes, -> case 2
29 RCL 04 a2
30 X<=Y? a2 <= pi ?
31 GTO 01 yes, -> case1
32 X<>Y no, it's case 3
33 XEQ 02 first between [pi, a2]
34 STO O partial result
35 RCL 03 a1
36 PI now between [a1, pi]
37 XEQ 02
38 RCL O
39 +
40 RTN
41 GTO 00
42 LBL 01
43 RCL 04 a2
44 RCL 03 a1
45 LBL 02
46 XEQ 05
47 STO M f1
48 X<>Y a2
49 XEQ 05
50 STO N f2
51 RCL M
52 -
53 2
54 / (f2 - f1) /2
55 ENTER^
56 SIN
57 RCL M f1
58 RCL N f2
59 +
60 2
61 / (f2 + f1) /2
62 COS
63 *
64 RCL 00 e
65 *
66 -
67 RCL 01 a
68 X^2 a^2
69 *
70 1
71 RCL 00
72 X^2
73 - 1-e^2
74 SQRT
75 *
76 ABS
77 RTN
78 GTO 00
79 LBL 05
80 COS cos a
81 RCL 00
82 X<>Y
83 + e + cos a
84 LASTX cos a
85 RCL 00 e
86 * e.cos a
87 1
88 +
89 / cos f = (e + cos a)/(1 + e.cos a)
90 ACOS f
91 END
(*) the symbol "<)" is for the angle character.
Example: calculate the area swept between a1 = pi/4 and a2 = 3.pi/4, if the ellipse parameters are a= 2, b= 3
Solution: A = 1.989554087
Once the area is obtained, and knowing the period of the orbiting (T), it's straightforward to determine the time taken by the planet to travel between the two positions, with the direct application of Kepler's 2nd. law:
t = T . Area / pi.a.b
Cheers,
ÁM