Post Reply 
Ellipsoid surface area
05-28-2021, 12:37 AM
Post: #1
Ellipsoid surface area
semi-axis, a ≥ b ≥ c
cos(φ) = c/a
m = (a²(b²-c²)) / (b²(a²-c²))
I = cos(φ)^2 * elliptf(φ,m) + sin(φ)^2 * ellipte(φ,m)

→ S = 2*pi*(c² + a*b/sin(φ)*I)

HP Prime does not have elliptf(ϕ,m) and ellipte(ϕ,m), we might as well combine I, as 1 integral

Let t = sin(ϕ)^2 = 1-(c/a)^2, s = √t = sin(ϕ):

I = (1-t) * ∫(1/√((1-x²)*(1-m*x²)), x=0..s) + t * ∫(√((1-m*x²)/(1-x²)), x=0..s)
I = ∫((1-m*t*x²) / √((1-x²)*(1-m*x²)), x=0..s)

Code:
#cas
ellipsoid_area(a,b,c)
BEGIN
    local m, t, s, x;
    c,b,a := sort(a,b,c);   // a >= b >= c
    t := 1.-(c/a)^2;
    if t==0 then return 4.*pi*a*c END
    s := sqrt(t);           // sin(acos(c/a))
    m := (b+c)*(b-c) / (b*b*t);
    m := int((1-m*t*x*x)/sqrt((1-x*x)*(1-m*x*x)), x, 0., s);
    return 2.*pi*(c*c + a*b/s*m);
END
#end

CAS> ellipsoid_area(1.1,1.2,4.7)             → 54.6901240998

Example from The Surface Area Of An Ellipsoid, A. Dieckmann, Universität Bonn, July 2003

A good estimate with Thomsen formula, rel. err ≤ 1.061%

CAS> ellipsoid_area_est(a,b,c) := 4*pi*mean(([a*b,a*c,b*c]).^p)^(1/p) | p = 1.6075

CAS> ellipsoid_area_est(1.1,1.2,4.7)       → 54.4952199737
Find all posts by this user
Quote this message in a reply
05-28-2021, 08:05 PM
Post: #2
RE: Ellipsoid surface area
Sorting the semi-axis may not be necessary.
This version removes semi-axis sorting step.
Also, test for perfect sphere is also eliminated.

Technically, with ellipsoid symmetry, input can be in any order.
However, we wanted to avoid catastrophic cancellation of terms.

As long as user input is in sorted order (or reverse sorted), it will give good result.
In other words, I am basing length relative to b.

Code:
#cas
ellipsoid_area(a,b,c)
BEGIN
    local t1, t2, x;
    t1 := 1.-(b/a)^2;
    t2 := 1.-(b/c)^2;
    t2 := int((1-t1*t2*x*x)/sqrt((1-t1*x*x)*(1-t2*x*x)), x, 0., 1);
    return 2.*pi*(b*b+t2*a*c);
END
#end

Note the t1, t2 symmetry. Swapping (a,c) returns same result.

CAS> ellipsoid_area(1,1,1)       → 12.5663706144 (4*pi)
CAS> ellipsoid_area(1,2,3)       → 48.8821463026
CAS> ellipsoid_area(1,3,2)       → 48.8821463026
CAS> ellipsoid_area(3,1,2)       → 48.8821463026

I confirmed numbers with keisan Online Calculator
For a,b,c = 3,2,1: S = 48.8821463025820596957
Find all posts by this user
Quote this message in a reply
05-30-2021, 06:04 PM (This post was last modified: 08-05-2022 04:57 PM by Albert Chan.)
Post: #3
RE: Ellipsoid surface area
This version use Carlson symmetric elliptic integrals, to avoid integration.

\(F\!\left(\phi, m\right) = \sin(\phi) R_F\!\left(\cos^{2}\!\left(\phi\right), 1 - m \sin^{2}\!\left(\phi\right), 1\right)\)

\(F\!\left(\phi, m\right) - E\!\left(\phi, m\right) = \frac{1}{3} m \sin^{3}\!\left(\phi\right) R_D\!\left(\cos^{2}\!\left(\phi\right), 1 - m \sin^{2}\!\left(\phi\right), 1\right)\)

We force R arguments to get close, with Duplication theorems.
Assume x,y,z > 0 (x or y can be zero, but not both)
Let λ = √(xy) + √(yz) + √(xz) > 0:

RF(x,y,z) = RF(x+λ, y+λ, z+λ)*2
RD(x,y,z) = RD(x+λ, y+λ, z+λ)*2 + 3/(√(z)*(z+λ))

If arguments close enough, we can stop.

RF(x,x,x) = x^(-1/2)
RD(x,x,x) = x^(-3/2)

From AM-GM inequality, we have (x+y)/2 ≥ √(xy)

λ = √(xy) + √(yz) + √(xz) ≤ (x+y)/2 + (y+z)/2 + (x+z)/2 = x+y+z

Code assumed λ*(1+1e-12) > x+y+z, we consider x,y,z close enough.

Code:
#cas
RFD(x,y,z)  // return RF(x,y,z), RD(x,y,z)
BEGIN
    LOCAL rx,ry,rz,k;
    rx := sqrt(approx(x));
    ry := sqrt(approx(y));
    rz := sqrt(approx(z));
    k := rx*(ry+rz) + ry*rz;
    rx := k*(1+1e-12) > x+y+z; // x,y,z close ?
    IF rx THEN rx:=sqrt(3/k); return rx, 3*rx/k END 
    rx, ry := RFD(x+k, y+k, z+k);
    return rx*2, ry*2 + 3/(rz*(z+k));
END;

ellipsoid_area(a,b,c)
BEGIN
    LOCAL rf, rd, k;
    rf := approx(b/a)^2;
    rd := approx(b/c)^2;
    k := (1-rf)*(1-rd)/3;
    rf, rd := RFD(rf,rd,1);
    return 2.*pi*(b*b+(rf-k*rd)*a*c);
END;    
#end

>>> mp.pretty = True
>>> elliprf(1,2,3), elliprd(1,2,3) # see https://mpmath.org/doc/current/functions/elliptic.html
(0.726945935468908, 0.290460281028991)

CAS> RFD(1,2,3)                     → [0.726945935469, 0.290460281029]
CAS> ellipsoid_area(1,2,3)       → 48.8821463026

Numerical Recipes, Elliptic Integrals: https://core.ac.uk/download/pdf/211378009.pdf
HP-41 Module: Elliptic Functions and Orthogonal Polynomials, ellipsoid area example (page 9), using RG
Find all posts by this user
Quote this message in a reply
05-31-2021, 01:38 AM
Post: #4
RE: Ellipsoid surface area
(05-30-2021 06:04 PM)Albert Chan Wrote:  HP-41 Module: Elliptic Functions and Orthogonal Polynomials, ellipsoid area example (page 9), using RG

From formula, S = 4*pi*RG((a*b)^2, (a*c)^2, (b*c)^2), we reverse it, for definition of RG:
Code:
#cas
RG(x,y,z)
BEGIN
    local ab, ac, bc;
    ab, ac, bc := map(sqrt, [x,y,z]);
    x, z := RFD(x,z,y); // relative to b^2
    z *= (bc-ac)*(bc+ac)*(ab-ac)*(ab+ac)/3;
    return (ab*bc/ac + (x*y-z))/2;
END;  
#end

CAS> ellipsoid_area(2,4,9)               → 283.427384268
CAS> 4*π*RG(8^2,18^2,36^2)       → 283.427384268
Find all posts by this user
Quote this message in a reply
08-05-2022, 04:40 PM (This post was last modified: 08-06-2022 09:23 PM by Albert Chan.)
Post: #5
RE: Ellipsoid surface area
Extended RFD(x,y,z) for complex numbers, replaced comparison test with equality (with fuzzy factor)

Code:
RFD(x,y,z)  // return RF(x,y,z), RD(x,y,z)
BEGIN
    LOCAL rx,ry,rz,k;
    rx := sqrt(approx(x));
    ry := sqrt(approx(y));
    rz := sqrt(approx(z));
    k := rx*(ry+rz) + ry*rz;
    rx := (x+y+z-k)*.01; // eps
    ry := (rx+k-k == 0); // x,y,z close?
//  print(k, k+20*rx, k+100*rx);
    IF ry THEN k:=sqrt(3/(k+20*rx)); return [k, k^3] END 
    rx, ry := RFD((x+k)/4, (y+k)/4, (z+k)/4);
    return [rx, ry/4 + 3/(rz*(z+k))];
END;

Previously, recursed arguments quadrupled in size.
Arguments (x,y,z) now scaled to be similar in size.

If we uncomment print(...) code, we get (x,y,z), convergence trend.
If x==y==z, then (k = √(xy)+√(yz)+√(xz)) == (x+y+z)

CAS> ellipsoid_area(2,4,9)      → 283.427384268
Code:
k              k+(x+y+z-k)/5  (x+y+z)
3.33333333333  3.70617283951  5.1975308642
3.65931727734  3.68733036508  3.79938271605
3.68498213244  3.68685243336  3.69433363702
3.68672507777  3.68684406393  3.68732000859
3.68683645911  3.68684392938  3.68687381047
3.68684345985  3.68684392727  3.68684579695
3.68684389801  3.68684392723  3.68684404412
3.68684392541  3.68684392723  3.68684393454
3.68684392712  3.68684392723  3.68684392769
3.68684392723  3.68684392723  3.68684392726
3.68684392723  3.68684392723  3.68684392723
3.68684392723  3.68684392723  3.68684392723

k seems to converge 4x faster, than (x+y+z), but from opposite side.
Extrapolate this, we have (k + (x+y+z-k)/5) converge at doubled rate.

Trivia: HP Prime float equality is really approximately equal (about 12 digits)

CAS> 9*(1+1e-12)            → 9.00000000001
CAS> Ans>9., Ans==9.      → [1, 1]

If our convergence test were "rx+k == k", we may lose needed recursions.
Test is written as "rx+k-k == 0", to ensure what we meant, bits equality.
Find all posts by this user
Quote this message in a reply
08-05-2022, 04:51 PM
Post: #6
RE: Ellipsoid surface area
Example: Elliptic Integral of the First Kind (note: we use elliptic parameter m=k^2)

>>> from mpmath import *
>>> ellipf(2,3)
1.0010773804561062-1.490278044744527j

RFD(x,y,z) can do both (F, E) at the same time.
Code:
FE(t,m) := // return [F(t,m), E(t,m)-F(t,m)]
BEGIN
t := sin(t);
m *= t*t;
return RFD(1-t*t,1-m,1) .* [t, -m*t/3];
END;

Formula is good for angle = 0 .. pi/2

CAS> k := FE(pi/2,3)[1]      → 1.00107738046-1.17142008415*i

// 2 = 2*pi/2 + (2-pi)
CAS> 2*k + FE(2-pi,3)[1]                             → 1.00107738046-1.49027804474*i

// 2 = pi/2 + (2-pi/2)
CAS> k + FE(2-pi/2,3/(3-1))[1]/sqrt(1-3)      → 1.00107738046-1.49027804474*i
Find all posts by this user
Quote this message in a reply
Post Reply 




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