Ellipsoid surface area - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: HP Prime Software Library (/forum-15.html) +--- Thread: Ellipsoid surface area (/thread-17013.html) Ellipsoid surface area - Albert Chan - 05-28-2021 12:37 AM 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 RE: Ellipsoid surface area - Albert Chan - 05-28-2021 08:05 PM 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 RE: Ellipsoid surface area - Albert Chan - 05-30-2021 06:04 PM 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 RE: Ellipsoid surface area - Albert Chan - 05-31-2021 01:38 AM (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 RE: Ellipsoid surface area - Albert Chan - 08-05-2022 04:40 PM 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. RE: Ellipsoid surface area - Albert Chan - 08-05-2022 04:51 PM 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