282 lines
8.1 KiB
Plaintext
282 lines
8.1 KiB
Plaintext
/*
|
|
Implementation of methods given in
|
|
B. C. Carlson
|
|
Computation of elliptic integrals
|
|
Numerical Algorithms 10, 13-26 (1995)
|
|
|
|
Copyright (c) Charles Karney (2009-2013) <charles@karney.com> and
|
|
licensed under the MIT/X11 License. For more information, see
|
|
https://geographiclib.sourceforge.io/
|
|
*/
|
|
|
|
/* fpprec:120$ Should be set outside */
|
|
etol:0.1b0^fpprec$ /* For Carlson */
|
|
tolRF : (3 * etol)^(1/8b0)$
|
|
tolRD : (etol/5)^(1/8b0)$
|
|
tolRG0 : 2.7b0 * sqrt(etol)$
|
|
tolJAC:sqrt(etol)$ /* For Bulirsch */
|
|
pi:bfloat(%pi)$
|
|
atan2x(y,x) := -atan2(-y, x)$ /* fix atan2(0b0,-1b0) = -pi */
|
|
rf(x,y,z) := block([a0:(x+y+z)/3, q,x0:x,y0:y,z0:z,
|
|
an,ln,xx,yy,zz,e2,e3,mul:1b0],
|
|
q:max(abs(a0-x),abs(a0-y),abs(a0-z))/tolRF,
|
|
an:a0,
|
|
while q >= mul * abs(an) do (
|
|
ln:sqrt(x0)*sqrt(y0)+sqrt(y0)*sqrt(z0)+sqrt(z0)*sqrt(x0),
|
|
an:(an+ln)/4,
|
|
x0:(x0+ln)/4,
|
|
y0:(y0+ln)/4,
|
|
z0:(z0+ln)/4,
|
|
mul:mul*4),
|
|
xx:(a0-x)/(mul*an),
|
|
yy:(a0-y)/(mul*an),
|
|
zz:-xx-yy,
|
|
e2:xx*yy-zz^2,
|
|
e3:xx*yy*zz,
|
|
(e3 * (6930 * e3 + e2 * (15015 * e2 - 16380) + 17160) +
|
|
e2 * ((10010 - 5775 * e2) * e2 - 24024) + 240240) /
|
|
(240240 * sqrt(an)))$
|
|
rd(x,y,z) := block([a0:(x+y+3*z)/5, q,x0:x,y0:y,z0:z,
|
|
an,ln,xx,yy,zz,e2,e3,e4,e5,s:0b0,mul:1b0],
|
|
q:max(abs(a0-x),abs(a0-y),abs(a0-z))/tolRD,
|
|
an:a0,
|
|
while q >= mul*abs(an) do (
|
|
ln:sqrt(x0)*sqrt(y0)+sqrt(y0)*sqrt(z0)+sqrt(z0)*sqrt(x0),
|
|
s:s+1/(mul*sqrt(z0)*(z0+ln)),
|
|
an:(an+ln)/4,
|
|
x0:(x0+ln)/4,
|
|
y0:(y0+ln)/4,
|
|
z0:(z0+ln)/4,
|
|
mul:4*mul),
|
|
xx:(a0-x)/(mul*an),
|
|
yy:(a0-y)/(mul*an),
|
|
zz:-(xx+yy)/3,
|
|
e2:xx*yy-6*zz^2,
|
|
e3:(3*xx*yy-8*zz^2)*zz,
|
|
e4:3*(xx*yy-zz^2)*zz^2,
|
|
e5:xx*yy*zz^3,
|
|
((471240 - 540540 * e2) * e5 +
|
|
(612612 * e2 - 540540 * e3 - 556920) * e4 +
|
|
e3 * (306306 * e3 + e2 * (675675 * e2 - 706860) + 680680) +
|
|
e2 * ((417690 - 255255 * e2) * e2 - 875160) + 4084080) /
|
|
(4084080 * mul * an * sqrt(an)) + 3 * s)$
|
|
|
|
/* R_G(x,y,0) */
|
|
rg0(x,y) := block([x0:sqrt(max(x,y)),y0:sqrt(min(x,y)),xn,yn,
|
|
t,s:0b0,mul:0.25b0],
|
|
xn:x0,
|
|
yn:y0,
|
|
while abs(xn-yn) > tolRG0 * xn do (
|
|
t:(xn+yn)/2,
|
|
yn:sqrt(xn*yn),
|
|
xn:t,
|
|
mul : 2*mul,
|
|
s:s+mul*(xn-yn)^2),
|
|
((x0+y0)^2/4 - s)*pi/(2*(xn+yn)) )$
|
|
|
|
rg(x, y, z) := (
|
|
if z = 0b0 then (z:y, y:0b0),
|
|
(z * rf(x, y, z) - (x-z) * (y-z) * rd(x, y, z) / 3
|
|
+ sqrt(x * y / z)) / 2)$
|
|
|
|
rj(x, y, z, p) := block([A0:(x + y + z + 2*p)/5,An,delta:(p-x) * (p-y) * (p-z),
|
|
Q,x0:x,y0:y,z0:z,p0:p,mul:1b0,mul3:1b0,s:0b0,X,Y,Z,P,E2,E3,E4,E5],
|
|
An : A0,
|
|
Q : max(abs(A0-x), abs(A0-y), abs(A0-z), abs(A0-p)) / tolRD,
|
|
while Q >= mul * abs(An) do block([lam,d0,e0],
|
|
lam : sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0),
|
|
d0 : (sqrt(p0)+sqrt(x0)) * (sqrt(p0)+sqrt(y0)) * (sqrt(p0)+sqrt(z0)),
|
|
e0 : delta/(mul3 * d0^2),
|
|
s : s + rc(1b0, 1b0 + e0)/(mul * d0),
|
|
An : (An + lam)/4,
|
|
x0 : (x0 + lam)/4,
|
|
y0 : (y0 + lam)/4,
|
|
z0 : (z0 + lam)/4,
|
|
p0 : (p0 + lam)/4,
|
|
mul : 4*mul,
|
|
mul3 : 64*mul3),
|
|
X : (A0 - x) / (mul * An),
|
|
Y : (A0 - y) / (mul * An),
|
|
Z : (A0 - z) / (mul * An),
|
|
P : -(X + Y + Z) / 2,
|
|
E2 : X*Y + X*Z + Y*Z - 3*P*P,
|
|
E3 : X*Y*Z + 2*P * (E2 + 2*P*P),
|
|
E4 : (2*X*Y*Z + P * (E2 + 3*P*P)) * P,
|
|
E5 : X*Y*Z*P*P,
|
|
((471240 - 540540 * E2) * E5 +
|
|
(612612 * E2 - 540540 * E3 - 556920) * E4 +
|
|
E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
|
|
E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
|
|
(4084080 * mul * An * sqrt(An)) + 6 * s)$
|
|
|
|
rd(x, y, z) := block([A0:(x + y + 3*z)/5,An,Q,x0:x,y0:y,z0:z,
|
|
mul:1b0,s:0b0,X,Y,Z,E2,E3,E4,E5],
|
|
An : A0,
|
|
Q : max(abs(A0-x), abs(A0-y), abs(A0-z)) / tolRD,
|
|
while Q >= mul * abs(An) do block([lam],
|
|
lam : sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0),
|
|
s : s + 1/(mul * sqrt(z0) * (z0 + lam)),
|
|
An : (An + lam)/4,
|
|
x0 : (x0 + lam)/4,
|
|
y0 : (y0 + lam)/4,
|
|
z0 : (z0 + lam)/4,
|
|
mul : 4*mul),
|
|
X : (A0 - x) / (mul * An),
|
|
Y : (A0 - y) / (mul * An),
|
|
Z : -(X + Y) / 3,
|
|
E2 : X*Y - 6*Z*Z,
|
|
E3 : (3*X*Y - 8*Z*Z)*Z,
|
|
E4 : 3 * (X*Y - Z*Z) * Z*Z,
|
|
E5 : X*Y*Z*Z*Z,
|
|
((471240 - 540540 * E2) * E5 +
|
|
(612612 * E2 - 540540 * E3 - 556920) * E4 +
|
|
E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
|
|
E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
|
|
(4084080 * mul * An * sqrt(An)) + 3 * s)$
|
|
|
|
rc(x, y):=if x < y then atan(sqrt((y - x) / x)) / sqrt(y - x)
|
|
else ( if x = y and y > 0b0 then 1 / sqrt(y) else
|
|
atanh( if y > 0b0 then sqrt((x - y) / x)
|
|
else sqrt(x / (x - y)) ) / sqrt(x - y) )$
|
|
|
|
/* k^2 = m */
|
|
ec(k2):=2*rg0(1b0-k2,1b0)$
|
|
kc(k2):=rf(0b0,1b0-k2,1b0)$
|
|
dc(k2):=rd(0b0,1b0-k2,1b0)/3$
|
|
pic(k2,alpha2):=if alpha2 # 0b0 then
|
|
kc(k2)+alpha2*rj(0b0,1b0-k2,1b0,1b0-alpha2)/3
|
|
else kc(k2)$
|
|
gc(k2,alpha2):=if alpha2 # 0b0 then (if k2 # 1b0 then
|
|
kc(k2) + (alpha2-k2)*rj(0b0,1b0-k2,1b0,1b0-alpha2)/3
|
|
else rc(1b0,1b0-alpha2))
|
|
else ec(k2)$
|
|
hc(k2,alpha2):=if alpha2 # 0b0 then (if k2 # 1b0 then
|
|
kc(k2) - (1b0-alpha2)*rj(0b0,1b0-k2,1b0,1b0-alpha2)/3
|
|
else rc(1b0,1b0-alpha2))
|
|
else kc(k2) - dc(k2)$
|
|
|
|
/* Implementation of methods given in
|
|
|
|
Roland Bulirsch
|
|
Numerical Calculation of Elliptic Integrals and Elliptic Functions
|
|
Numericshe Mathematik 7, 78-90 (1965)
|
|
|
|
*/
|
|
|
|
sncndn(x,mc):=block([bo, a, b, c, d, l, sn, cn, dn, m, n],
|
|
local(m, n),
|
|
if mc # 0 then (
|
|
bo:is(mc < 0b0),
|
|
if bo then (
|
|
d:1-mc,
|
|
mc:-mc/d,
|
|
d:sqrt(d),
|
|
x:d*x),
|
|
dn:a:1,
|
|
for i:0 thru 12 do (
|
|
l:i,
|
|
m[i]:a,
|
|
n[i]:mc:sqrt(mc),
|
|
c:(a+mc)/2,
|
|
if abs(a-mc)<=tolJAC*a then return(false),
|
|
mc:a*mc,
|
|
a:c
|
|
),
|
|
x:c*x,
|
|
sn:sin(x),
|
|
cn:cos(x),
|
|
if sn#0b0 then (
|
|
a:cn/sn,
|
|
c:a*c,
|
|
for i:l step -1 thru 0 do (
|
|
b:m[i],
|
|
a:c*a,
|
|
c:dn*c,
|
|
dn:(n[i]+a)/(b+a),
|
|
a:c/b
|
|
),
|
|
a:1/sqrt(c*c+1b0),
|
|
sn:if sn<0b0 then -a else a,
|
|
cn:c*sn
|
|
),
|
|
if bo then (
|
|
a:dn,
|
|
dn:cn,
|
|
cn:a,
|
|
sn:sn/d
|
|
)
|
|
) else /* mc = 0 */ (
|
|
sn:tanh(x),
|
|
dn:cn:sech(x)
|
|
/* d:exp(x), a:1/d, b:a+d, cn:dn:2/b,
|
|
if x < 0.3b0 then (
|
|
d:x*x*x*x,
|
|
d:(d*(d*(d*(d+93024b0)+3047466240b0)+24135932620800b0)+
|
|
20274183401472000b0)/60822550204416000b0,
|
|
sn:cn*(x*x*x*d+sin(x))
|
|
) else
|
|
sn:(d-a)/b */
|
|
),
|
|
[sn,cn,dn]
|
|
)$
|
|
|
|
Delta(sn, k2):=sqrt(1 - k2 * sn*sn)$
|
|
|
|
/* Versions of incomplete functions in terms of Jacobi elliptic function
|
|
with u = am(phi) real and in [0,K(k2)] */
|
|
|
|
eirx(sn,cn,dn,k2,ec):=block([t,cn2:cn*cn,dn2:dn*dn,sn2:sn*sn,ei,kp2:1b0-k2],
|
|
ei : (if k2 <= 0b0 then
|
|
rf(cn2, dn2, 1b0) - k2 * sn2 * rd(cn2, dn2, 1b0) / 3 else
|
|
(if kp2 >= 0b0 then
|
|
kp2 * rf(cn2, dn2, 1b0) +
|
|
k2 * kp2 * sn2 * rd(cn2, 1b0, dn2) / 3 +
|
|
k2 * abs(cn) / dn else
|
|
- kp2 * sn2 * rd(dn2, 1b0, cn2) / 3 + dn / abs(cn) ) ),
|
|
ei : ei * abs(sn),
|
|
if cn < 0b0 then ei : 2 * ec - ei,
|
|
if sn < 0 then ei : -ei,
|
|
ei)$
|
|
|
|
dirx(sn, cn, dn, k2, dc):=block(
|
|
[di:abs(sn) * sn*sn * rd(cn*cn, dn*dn, 1b0) / 3],
|
|
if cn < 0b0 then di : 2 * dc - di,
|
|
if sn < 0b0 then di : -di,
|
|
di)$
|
|
|
|
hirx(sn, cn, dn, k2, alpha2, hc):=block(
|
|
[cn2 : cn*cn, dn2 : dn*dn, sn2 : sn*sn, hi, alphap2:1b0-alpha2],
|
|
hi : abs(sn) * (rf(cn2, dn2, 1b0) -
|
|
alphap2 * sn2 *
|
|
rj(cn2, dn2, 1b0, cn2 + alphap2 * sn2) / 3),
|
|
if cn < 0 then hi : 2 * hc - hi,
|
|
if sn < 0 then hi : -hi,
|
|
hi)$
|
|
|
|
deltae(sn, cn, dn, k2, ec):=(
|
|
if cn < 0 then (cn : -cn, sn : -sn),
|
|
eirx(sn, cn, dn, k2, ec) * (pi/2) / ec - atan2x(sn, cn))$
|
|
deltad(sn, cn, dn, k2, dc):=(
|
|
if cn < 0 then (cn : -cn, sn : -sn),
|
|
dirx(sn, cn, dn, k2, ec) * (pi/2) / dc - atan2x(sn, cn))$
|
|
deltah(sn, cn, dn, k2, alpha2, hc):=(
|
|
if cn < 0 then (cn : -cn, sn : -sn),
|
|
hirx(sn, cn, dn, k2, alpha2, hc) * (pi/2) / hc - atan2x(sn, cn))$
|
|
|
|
einv(x, k2, ec):=block(
|
|
[n : floor(x / (2 * ec) + 0.5b0), phi, eps : k2/(sqrt(1b0-k2) + 1b0)^2,
|
|
err:1b0],
|
|
x : x - 2 * ec * n,
|
|
phi : pi * x / (2 * ec),
|
|
phi : phi - eps * sin(2 * phi) / 2,
|
|
while abs(err) > tolJAC do block([sn:sin(phi), cn:cos(phi), dn],
|
|
dn : Delta(sn, k2),
|
|
err : (eirx(sn, cn, dn, k2, ec) - x)/dn,
|
|
phi : phi - err),
|
|
n * pi + phi)$
|
|
|
|
deltaeinv(stau, ctau, k2, ec) := block([tau],
|
|
if ctau < 0 then (ctau : -ctau, stau : -stau),
|
|
tau : atan2x(stau, ctau),
|
|
einv(tau * ec / (pi/2), k2, ec) - tau)$
|