Post Reply 
Perimeter of Ellipse
01-23-2020, 01:40 PM (This post was last modified: 01-24-2020 05:04 PM by Albert Chan.)
Post: #21
RE: Perimeter of Ellipse
(01-21-2020 05:16 PM)Albert Chan Wrote:  \(\large p(a,b) = 2× \left( p({a+b \over 2}, \sqrt{ab}) - {\pi a b \over AGM(a,b)}\right)\)

In other words, calculate ellipse perimeter from a much less eccentric ellipse.

We could even do it in a novel way, from the *more* eccentric ellipse.
Below, LHS is the more eccentric ellipse, with \({a_0 + b_0 \over 2} = a, \sqrt{a_0 b_0} = b\)

\( p(a_0,b_0) = 2× \left( p(a,b) - {\pi b^2 \over AGM(a,b)}\right)\)

Using my Casio FX-115MS, AGM2 method for Halley's comet, numbers from here

A = 2667950000
B = 678281900
C = A²
D = B²
E = 0.5

F=B-A : C=C-EF² : B=√AB : A=A+F/2 : E=2E

Quote:\(p(a,b) = {p(a_0,b_0) \over 2} + {\pi b^2 \over AGM(a,b)}\)

Keep pressing "=" until F=0 (so that A,B,C all converged)

\(\large\pi\)(C+D)/B       → ellipse perimeter = 11,464,318,984.1 km

This was similar to my lua code, with C=(A²+B²)/2, E=0.25, so that converged C/B = "radius".
Find all posts by this user
Quote this message in a reply
06-05-2020, 03:28 AM (This post was last modified: 06-06-2020 11:33 AM by Albert Chan.)
Post: #22
RE: Perimeter of Ellipse
From previous 2 posts, we can get some nice Elliptic Integral Identity:
Below, all Elliptic functions uses modulus argument.

Eccentricity of ellipse, e^2 = 1 - (b/a)^2
We proved "next" ellipse (a' = (a+b)/2 , b' = sqrt(a*b)), has eccentrity of h = (a-b)/(a+b)

K(e) = pi/2 / agm(1,b/a) = pi/2 / agm(1+e,1-e)
K(h) = pi/2 / agm(1+h,1-h) = pi/2 / agm(1,b/a) / (1+h)

K(e) = (1+h) K(h)

(01-21-2020 05:16 PM)Albert Chan Wrote:  \(\large p(a,b) = 2× \left( p({a+b \over 2}, \sqrt{ab}) - {\pi a b \over AGM(a,b)}\right)\)

With ellipse_perimeter = 4 a E(e), we have:

4 a E(e) = 4 (a+b) E(h) - 4 b K(e) = 4 (a+b) E(h) - 4 b (1+h) K(h)

E(e) = 2/(h+1) * E(h) - (1-h) * K(h)
Find all posts by this user
Quote this message in a reply
06-05-2020, 08:09 PM (This post was last modified: 06-06-2020 05:14 PM by Albert Chan.)
Post: #23
RE: Perimeter of Ellipse
We can get E(e),K(e) by recursively going for E(h),K(h), like this.
Note: code use parameter m, not modulus k, m = k^2
Note: code quit when m is small enough for taylor linear term estimate.

E(m=ε) ≈ pi/2 - (pi/8)*ε ,     K(m=ε) ≈ pi/2 + (pi/8)*ε

Code:
from cmath import sqrt, pi
    
def EK(m, verbal=False, c=pi/2):
    e = abs(m)
    if 1+e*e == 1: m *= 0.25*c; return c-m, c+m
    b = sqrt(1-m)
    if not b: return 1, 99  # K(1) assumed 99 (should be inf)
    h = m / ((b+2)*b+1)     # = (1-b)/(1+b)
    e, k = EK(h*h, verbal)
    hk = h*k
    e = e-k + hk + b*e      # = 2/(1+h)*e - (1-h)*k
    k = hk + k              # = (1+h)*k
    if verbal: print 'm = %s\n\tE(m) = %s\n\tK(m) = %s' % (m, e, k)
    return e, k

>>> e, k = EK(0.96, verbal=True)
m = (2.89332317725e-05+0j)
      E(m) = (1.57078496468+0j)
      K(m) = (1.57080768903+0j)
m = (0.0212862362522+0j)
      E(m) = (1.56240357945+0j)
      K(m) = (1.57925700384+0j)
m = (0.444444444444+0j)
      E(m) = (1.3781039379+0j)
      K(m) = (1.80966749549+0j)
m = 0.96
      E(m) = (1.05050222698+0j)
      K(m) = (3.01611249248+0j)
>>> 4 * 50 * abs(e)                       # ellipse_perimeter(10,50), error = -3 ULP
210.10044539689011

>>> e, k = EK(2, verbal=True)        # see https://www.hpmuseum.org/forum/thread-15...#pid132745
m = (5.57959210499e-05-0j)
      E(m) = (1.57077441556+0j)
      K(m) = (1.57081823849+0j)
m = (0.0294372515229-0j)
      E(m) = (1.55917174457+0j)
      K(m) = (1.58255172722+0j)
m = (-1-0j)
      E(m) = (1.91009889451+0j)
      K(m) = (1.31102877715+0j)
m = 2
      E(m) = (0.599070117368+0.599070117368j)
      K(m) = (1.31102877715-1.31102877715j)
Find all posts by this user
Quote this message in a reply
06-06-2020, 05:12 PM (This post was last modified: 06-09-2020 12:55 AM by Albert Chan.)
Post: #24
RE: Perimeter of Ellipse
Improvement to my previous EK(m)
When m is small, the term E(m) - K(m) is losing significant digits.
Worse, catastrophic cancellation occurs at the base of recursion.

EKmc(m) returns E(m)-c, K(m)-c, where c=pi/2, similarly to expm1(x) = exp(x) - 1
EK(m) is now a simple wrapper for EKmc(m), adding back the c's.

Code:
from cmath import sqrt, pi
c = pi/2

def EKmc(m):
    'return E(m)-c, K(m)-c, where c=pi/2, m!=1'
    if 1+abs(m) == 1: m *= 0.25*c; return -m, m
    b = sqrt(1-m)
    h = m / ((b+2)*b+1)     # = (1-b)/(1+b)
    e, k = EKmc(h*h)
    hk, hc = h*k, h*c
    return (hk-k) + (e+b*e) - b*hc, (hk+k) + hc

def EK(m):
    'return E(m), K(m), assumed K(1)=99 (should be inf)'
    if abs(m-1)==0: return 1, 99
    e, k = EKmc(m)
    return e+c, k+c

>>> e, k = EK(0.96)
>>> 4 * 50 * abs(e)                       # ellipse_perimeter(10,50), error = 0 ULP
210.10044539689002
Find all posts by this user
Quote this message in a reply
08-01-2020, 12:31 PM
Post: #25
RE: Perimeter of Ellipse
(01-21-2020 05:16 PM)Albert Chan Wrote:  
(01-19-2020 11:00 PM)Albert Chan Wrote:  \(\large p(a,b) = 2× \left( p({a+b \over 2}, \sqrt{ab}) - {\pi a b \over AGM(a,b)}\right)\)

In other words, calculate ellipse perimeter from a much less eccentric ellipse.

Amazingly, the less eccentric ellipse have the eccentricity of h
[Image: ellipaxes.gif]
If we are at focus F2, and we let distance a0 = AF2, b0 = BF2, eccentricity of ellipse is also h0 Smile

ellipse major-axis a = OB = (a0+b0)/2
ellipse minor-axis b = OC = √((CF2)² - (OF2)²) = √(a² - (a-b0)²) = √((2a-b0)*b0) = √(a0*b0)

This matched ellipse perimeter AGM formula !

Example, lets try the earth-moon orbit, average orbital distance 384748 km

a0 = 406731 km
b0 = 364397 km

e = h0 = (a0-b0)/(a0+b0) = 42334/771128 ≈ 0.0549
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)