# 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
