# The complete elliptic integral of the first kind in terms of the modulus "k" . # Reference: Abramowitz and Stegun, chapter 17 # Method: Arithmetic-Geometric mean # G. Fee Nov. 1990 `evalf/EllipticK` := proc(kk) local a,b,c,k,e,g,t; k := evalf(kk); if not type(k,numeric) or k < 0 or k > 1 then RETURN('EllipticK'(k)) fi; b := 1-k; if b=0 then ERROR(`singularity encountered`) else e := Float(1,-Digits); g := Digits; Digits := Digits + length(Digits) + 1; a := 1.; c := k; Digits := 2*Digits; b := 1-k^2; Digits := iquo(Digits,2); b := b^(1/2); while c>e do t := a-b; c := .5*t; t := a*b; a := b+c; b := t^(1/2); od; t := evalf(.5*Pi/a); Digits := g; evalf(t); fi end: # The complete elliptic integral of the second kind in terms of the # modulus "k" . # Reference: Abramowitz and Stegun, chapter 17 # Method: Arithmetic-Geometric mean # G. Fee Nov. 1990 `evalf/EllipticE` := proc(kk) local t,a,b,c,k,e,p,s,g; k := evalf(kk); if not type(k,numeric) or k < 0 or k > 1 then RETURN('EllipticE'(k)) fi; b := 1-k; if b=0 then 1. else e := Float(1,-Digits); g := Digits; Digits := Digits + length(Digits) + 1; p := 1; a := 1.; c := k; s := k^2; Digits := 2*Digits; b := 1-k^2; Digits := iquo(Digits,2); b := b^(1/2); while c>e do t := a-b; c := .5*t; t := a*b; a := b+c; b := t^(1/2); p := 2*p; s := s + p*c^2; od; t := evalf(.5*Pi/a*(1-.5*s)); Digits := g; evalf(t); fi end: # save `evalf/EllipticK`, `EllipticK.m`; # save `evalf/EllipticE`, `EllipticE.m`; # quit