Kepler's 2nd. Law

[Ángel Martin]

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.

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

[toml_12953]

(02-12-2018 09:53 AM)Ángel Martin Wrote:  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.

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

For the parameters a=2, b=3 and a1=.785398..., a2=2.35619...

(Did I get the parameters right?)

I get 5.89676233948396 which agrees with the calculator at

http://keisan.casio.com/exec/system/1343722259

Code:

DECLARE EXTERNAL FUNCTION F
INPUT PROMPT "Enter semimajor axis a: ":a
INPUT PROMPT "Enter semiminor axis b: ":b
INPUT PROMPT "Enter Theta0: ":theta0
INPUT PROMPT "Enter Theta1: ":theta1
LET S = F(a,b,Theta1) - F(a,b,Theta0)
PRINT "Area =";S
END
EXTERNAL FUNCTION F(a,b,t)
LET F = a*b/2*(t-ATN((b-a)*SIN(2*t)/(b+a+(b-a)*COS(2*t))))
END FUNCTION

[SlideRule]

Angel

The 143 page publication Introduction to Orbital Flight Planning (I) (NASA-CE-165052) may be of interest as it includes a well-documented {equations, diagrams, etc.}discussion with respect to orbital mechanics & Kepler as well as an HP-67 based program for same. The discussion starts with elliptic orbits on page 51 then transitions to Kepler on page 54 with an HP-67 program on page 57.

BEST!
SlideRule

[Ángel Martin]

(02-13-2018 01:07 AM)toml_12953 Wrote:  

(02-12-2018 09:53 AM)Ángel Martin Wrote:  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

For the parameters a=2, b=3 and a1=.785398..., a2=2.35619...

(Did I get the parameters right?)

I get 5.89676233948396 which agrees with the calculator at

http://keisan.casio.com/exec/system/1343722259

Mmm... something's fishy. I consistently get 1.989554087 using a1 = pi/4 (45 degrees) and a2 = 3.pi/4 (135 degrees).
Besides, it's not possible to get the same as the web applet - which uses centered sectors, not focus-centered ones.

In fact, using the angles above the applet shows an Area result of 3.5280156212854

Could you double check?

[Ángel Martin]

(02-13-2018 02:01 AM)SlideRule Wrote:  The 143 page publication Introduction to Orbital Flight Planning (I) (NASA-CE-165052) may be of interest as it includes a well-documented {equations, diagrams, etc.}discussion with respect to orbital mechanics & Kepler as well as an HP-67 based program for same. The discussion starts with elliptic orbits on page 51 then transitions to Kepler on page 54 with an HP-67 program on page 57.

Great reference, thanks for the link. I'll check it out ASAP.

[Luigi Vampa]

(02-13-2018 02:01 AM)SlideRule Wrote:  Angel
The 143 page publication Introduction to Orbital Flight Planning (I) (NASA-CE-165052) may be of interest as it includes a well-documented {equations, diagrams, etc.}discussion with respect to orbital mechanics & Kepler as well as an HP-67 based program for same. The discussion starts with elliptic orbits on page 51 then transitions to Kepler on page 54 with an HP-67 program on page 57.
BEST!
SlideRule

Big thanks for the information!

[Ángel Martin]

I've prepared a set of functions and FOCAL drivers for an "Elliptic Applications" ROM, which also includes a basic "Orbital Mechanics for Dummies" section - quite an slippery subject, if I may say. Got the geometry and satellite positioning taken care of, but the velocities are giving some inconsistent result, still needs some work... Stay tuned!

[Ángel Martin]

The module is almost ready, here's the function QRG to whet your appetite:

Code:

XROM    Function     Description                      Input / Output
16,00   -ELLIP APPS  Section Header                   n/a
16,01   BRHM         Area of Cyclic Quadrilateral     Data in R01-R04
16,02   CIRCLE       Circle through three points      P1, P2, P3 in Registers R01-R06
16,03   "ECS+"       Driver for ECS                   Prompts for values
16,04   ECS          Elliptical Central Sector        {a,b,e} in R00-R02,  in X
16,05   EFS-         Ellliptical Left Focal Sector    {a,b,e} in R00-R02,  in R03-04
16,06   “EFS+“       Elliptical Right Focal Sector    {a,b,e} in R00-R02,  in R03-04
16,07   HERON        Area of Triangle from sides      Data in X,Y,Z
16,08   "F>C“        Focal to Central angles          {a,b,e} in R00-R02,  in X
16,09   "K2-“        Kepler 2nd. Law – Left Sector    Prompts for values
16,10   "K2+“        Kepler 2nd. Law – Right Sector   Prompts for values
16,11   "LC0"        Central sector Arc Length        Prompts for values
16,12   "LF+"        Right Focal sectorArc Length     Prompts for values
16,13   RAMA2        Ramanujan’s 2nd. Approx.         a,b, in Y,X
16,14   ZELIP1       Complex Elliptic Intg. 1st kind  z in (Z,Y), m in X
16,15   ZELIP2       Complex Elliptic Intg. 2nd kind  z in (Z,Y), m in X
16,16   -ORBITS 101  Section Header                   n/a
16,17   AZMTH        Azimuth angle                    {a,b,e} in R00-R02,  in X
16,18   E>M          Eccentric to Mean anomaly        (e, E) in Y,X
16,19   E>T          Eccentric to True Anomaly        (e, E) in Y, X
16,20   M>E          Mean to Eccentric anomaly        (e, M) in Y,X
16,21   T-E          True to Eccentric Anomaly        (e, T) in Y,X
16,22  “ORBIT“       Delta-V Orbit Simulator          Under program control
16.23   “PRJTL“      Projectile Trajectory            Under Program control
16.24   “TA<T>“      True Anomaly via Kepler          New Position w/Kepler equation
16.25   “TA+“        True Anomaly Direct Mode (+)     New Position using Right sectors
16.26   “TA-“        True Anomaly Direct Mode (-)     New Position using Left sectors
16.27   “<)T“        Subroutine for Solve             Under program control
16.28   “V<R>“       Velocity “Vis Vivas“             Velocity at distance – Left focus
16.29   “VR+“        Velocity at Position R+          Velocity at distance - right focus
16.30   “VR-“        Velocity at Position R-          Velocity at distance - Left focus
16.31   -AUX FNS     Section Header                   n/a
16.32   <)C0+        Central angle for Focal 90°      (a b,e) in R00-R02. Lifts stack
16.33   <)C-R        Central angle to Radius          {a,b,e} in R00-R02,  in X
16.34   <)C-XY       Central angle to coordinates     {a,b,e} in R00-R02,  in X
16.35   <)F0+        Right focal angle at x=0         (a,b,e) in R00-R02. Lifts stack
16.36   <)F0-        Left Focal angle at x=0          (a,b,e) in R00-R02. Lifts stack
16.37   <)F>R+       Right focal angle to radius      {a,b,e} in R00-R02,  in X
16.38   <)F>R-       Left Focal angle to radius       {a,b,e} in R00-R02,  in X
16.39   “?ab“        Promptsfor semi-axis             New value or R/S for current
16,40   “?<)-“       Prompts for left-angles          New value or R/S for current
16,41   “?<)+“       Prompts for right-angles         New value or R/S for current
16,42   “?P“         Prompts for Period               New value or R/S for current
16,43   1/R^2        Inverse of squared radius        {a,b,e} in R00-R02,  in X
16,44   QROOT        Quadratic Eq. Roots              c2, c1, c0 in Z,Y,X
15,45   QROUT        Shows roots                      Roots in X,Y,Z
15.46   THV          Tetrahedron Volume               Data in R00-R06
15.47   V(X)         Velocity at position X           Velocity at abccisa - Central
15.48   ZOUT         Shows complex value              Re,Im in X, Y
15.49   -SOLAR SYS   Section Header                   n/a
15.50   EARTH        Planet orbital data              Puts data in R00-R02 & R06
15.51   JUPITER      Planet orbital data              Puts data in R00-R02 & R06
15.52   MARS         Planet orbital data              Puts data in R00-R02 & R06
15.53   MERCURY      Planet orbital data              Puts data in R00-R02 & R06
15.54   MOON         Satellite orbital data           Puts data in R00-R02 & R06
15.55   NEPTUNE      Planet orbital data              Puts data in R00-R02 & R06
15.56   PLUTO        Planet orbital data              Puts data in R00-R02 & R06
15.57   SATURN       Planet orbital data              Puts data in R00-R02 & R06
15.58   URANUS       Planet orbital data              Puts data in R00-R02 & R06
15.59   VENUS        Planet orbital data              Puts data in R00-R02 & R06
15.60   +V1          Hohmann Transfer #1              T, R1, R2 in (Z,Y,X)
15,61   +V2          Hohmann Transfer #2              T, R1, R2 in (Z,Y,X)
15,62   1/2         Gravitational Parameter          Period in Y, semi-major in X
16,63   T()         Period from Gravitational data   Mass in Y, semi-major axis in X

[twdeckard]

Fantastic. I’ll look forward to the ROM.

I’ve been working on a package to determine real time altitude and azimuth from a satellites two-line-elements. I am hopelessly mired in the coordinate transformations.

I posted a copy of a ham radio article from ORBIT that does the same, but it wasn’t as enamored with the advantage module as I am.

It is a fun exercise to point out visible objects like the ISS and the Humanity Star and I have a working arduino version which points an antenna array.

https://youtu.be/VG5L5oP2uNk

As a conceit I have thought about plumbing an HP-IL/RS232 to it to place the ‘41 in charge. With the single precision math of the arduino the orbital elements age very quickly.

https://github.com/twdeckard/HP41CX-PLAN...PLAN41.rpn
(the only thing salvageable in the code would be the keps parser)

Thanks for creating and sharing this.

Todd

[SlideRule]

(03-11-2018 01:45 PM)twdeckard Wrote:  … I posted a copy of a ham radio article from ORBIT that does the same, but it wasn’t as enamored with the advantage module…
Todd

Was that article Tracking Satellites in Elliptical Orbits pgs 46-50, Ham Radio Magazine 1981 by Paul C. Bunnell?

Best!
SlideRule

[Ángel Martin]

(03-11-2018 01:45 PM)twdeckard Wrote:  I’ve been working on a package to determine real time altitude and azimuth from a satellites two-line-elements. I am hopelessly mired in the coordinate transformations.

I posted a copy of a ham radio article from ORBIT that does the same, but it wasn’t as enamored with the advantage module as I am.

It is a fun exercise to point out visible objects like the ISS and the Humanity Star and I have a working arduino version which points an antenna array.

https://youtu.be/VG5L5oP2uNk

As a conceit I have thought about plumbing an HP-IL/RS232 to it to place the ‘41 in charge. With the single precision math of the arduino the orbital elements age very quickly.

https://github.com/twdeckard/HP41CX-PLAN...PLAN41.rpn
(the only thing salvageable in the code would be the keps parser)

that sure looks like serious skills at work, very impressing!

I hope you don"t set your expectations too high, the ROM is quite basic in contents and doesn't delve deep into the subjects - more of a personal past-time than any other thing.

BTW would you rather have your program in a ROM format? I can transfer it to ROM if that makes any sense to you...

Best,
ÁM

[SlideRule]

(02-12-2018 09:53 AM)Ángel Martin Wrote:  … 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…
ÁM …

NOT better, but DG references nonetheless (for those with an interest) …

[attachment=5735] [attachment=5736] [attachment=5737] [attachment=5738] [attachment=5739]

… as well as these three articles from Surveying & Mapping magazine …
Reference Ellipsoids direct & indirect solution v40 (09-80)
Reference Ellipsoids comments on formulas v41 03'81
Reference Ellipsoids comments on article v41 09'81

… with a smidgen of applicability from the following articles …
Total Inverse solution for the geoDesic & great Elliptic
New Solutions to the Direct and Indirect Geodetic Problems on the Ellipsoid
Direct & Inverse Solutions of Geodetics of the Ellipsoid with Application of Nested Equations
Solving the Direct & Inverse Geodetic Problems on the Ellipsoid by Numerical Integration
Ellipsometer Data Analysis w a Small Programmable Desk Calculator

Ignore what doesn't facilitate the projected outcome, enjoy the rest.

BEST!
SlideRule

[twdeckard]

(03-11-2018 03:37 PM)Ángel Martin Wrote:  BTW would you rather have your program in a ROM format? I can transfer it to ROM if that makes any sense to you...

Best,
ÁM

Thank you for this kind offer! I have the nutstudio tools so if I ever get it to work I think I can spin it into a ROM of focal routines.

Take care
Todd

[SlideRule]

Tracking Satellites in Elliptic Orbits article {referenced above} attached:

[attachment=5740] [attachment=5741] [attachment=5742] [attachment=5743] [attachment=5744]

BEST!
SlideRule

[Ángel Martin]

Thanks for sharing this article, will study it closely!

Best,
ÁM

[Ángel Martin]

(03-11-2018 04:06 PM)twdeckard Wrote:  

(03-11-2018 03:37 PM)Ángel Martin Wrote:  BTW would you rather have your program in a ROM format? I can transfer it to ROM if that makes any sense to you...

Thank you for this kind offer! I have the nutstudio tools so if I ever get it to work I think I can spin it into a ROM of focal routines.

Cool. BTW looking at the code a couple of things caught my attention:

• This is definitely way more advanced and specific than the module I've put together.
• Are the contents of the ASCII files "ISS" and "SATREG" available? [
• What's MMDOT2 (appears in LBL "MA+ORB: ). It appears to be MCODE function...
  

Best,
ÁM

[twdeckard]

Greetings Ángel

Thank you for the kind words, remember, the program doesn't work yet. Candidly the best effort would be to scrap it and weld the TLE parser on the article that SlideRule posted.

The SATREG ASCII file is a list of variable mnemonics prefixed with a numerical type. Its a common pattern I used so that I can revisit an old program and have some hope of deciphering it. It is also why I squander precious registers on six character labels. I invariably have to scavenge those as I get to the limit of memory.

Here is the SATREG contents:

Code:

4^TMP
4^TMP
4^TMP
4^TMP        
4^SX
4^SY
4^SZ
4^OX
4^OY
4^OZ
4^VOx
4^VOy
4^VOz
4^RX
4^RY
4^RZ
4ELAPST
4CURORB
4OBJECT
1EPYEAR
2EPTIME
3INCLIN
2R.A.A.N
4ECCENT
2PERIGE
2MEANAN
3MMDOT1
3MMDOT2
4DECAYR
1ORBITN
6OBSLAT
6OBSLON
4LAMBDA
1SATRNG
4DOPPLR
6SATLAT
6SATLON
1AZIMUT
1ELEVAT
1OBSHGT
4MEANAT
4KD    
4KDP

The ISS ASCII file contains the two-line element (sometimes slang'd as a Kepler although it is not). They can be had from Celestrak or one of the other providers although NORAD is the source.

Code:

1 25544U 98067A   18071.70991576  .00002716  00000-0  48084-4 0  9995
2 25544  51.6420 139.6732 0001880 178.1916 325.8769 15.54223294103518

The MMDOT refers to the decay rate. There certainly isn't any MCODE in my vocabulary. I mixed in comments from the original BASIC program that inspired this effort. Maybe the comments looked like assembly. If it hurts my credibility, at one point, I was comparing output from a Commodore 64 to error check my '41.

Take care if hoping for some sensible output from the github repository. I had an ESD catastrophe during development and I am not sure it is even up to date. I only now have reconstituted all the keystrokes in the calculator and it does not even parse the TLE correctly. Most likely due to typos.

Even when working I only got as far as a correct ECF to ENU transformation. Everything else was wrong. Rather than simply "port" the BASIC I tried to go back to the actual coordinate transformations but it has gotten the better of me. I bought a "backup" HP41CX off of ebay and it had an Advantage module so I was looking for an excuse to do some matrix operations.

I will post something if I make any forward progress.

Best
Todd

[Ángel Martin]

(03-12-2018 10:23 PM)twdeckard Wrote:  Greetings Ángel

Thank you for the kind words, remember, the program doesn't work yet. Candidly the best effort would be to scrap it and weld the TLE parser on the article that SlideRule posted.
.......

Even when working I only got as far as a correct ECF to ENU transformation. Everything else was wrong. Rather than simply "port" the BASIC I tried to go back to the actual coordinate transformations but it has gotten the better of me. I bought a "backup" HP41CX off of ebay and it had an Advantage module so I was looking for an excuse to do some matrix operations.

I will post something if I make any forward progress.

Thanks for the clarification, it'll stop me from chasing a red herring ;-)

BTW you're not alone having used the C64 as an "advanced" parallel development tool. Many moons ago I did the same for the synchronous machine swing equation... aah those good times!

My module was completed but then I decided to add formulas for the velocity increments for coaxial elliptical orbit transfers... and I'm now polishing a couple of rough edges.

Best,
ÁM

[SlideRule]

(02-12-2018 09:53 AM)Ángel Martin Wrote:  … 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 …

Is this diagram …

[attachment=5756]

… representative?

SlideRule

[Ángel Martin]

(03-14-2018 05:54 PM)SlideRule Wrote:  

(02-12-2018 09:53 AM)Ángel Martin Wrote:  … 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 …

Is this diagram …
… representative?

Absolutely.