function [phix] = ellipticf(phi, k)
a0 = 1.0;
b0 = sqrt(1-k);
c0 = sqrt(k);
phi0 = phi;
ch = 1e-10;
i = 0;
while ch < abs(c0)
  kk = floor((phi0+pi/2)/pi);
  phi0 = phi0+atan(b0/a0*tan(phi0-kk*pi))+kk*pi;
  c0 = (a0-b0)/2;
  a0p = (a0+b0)/2;
  b0 = sqrt(a0*b0);
  a0 = a0p;
  i = i+1;
end
phix = phi0/(2^i)/a0;