function [sn,cn,dn]=egellip(h,phiin,m,ny)%,tol)
%Solve the elliptic integral from 0 to phi 

%if nargin<3, tol = eps; end
%if nargin<2, error('Not enough input arguments.'); end

a(1) = 1;
c(1) = sqrt(m);
b(1) = sqrt(1-m);

i = 1;
while abs(c(i)) > eps
    i = i + 1;
    a(i) = 0.5 * (a(i-1) + b(i-1));
    b(i) = sqrt(a(i-1)* b(i-1));
    c(i) = 0.5 * (a(i-1) - b(i-1));
end
n=i-1;
phiit(1)=phiin;
for j=1:n
    phijf=mod(phiit(j),pi);
    nte=(phiit(j)-phijf)/pi;
    tempf=atan((b(j)/a(j))*tan(phijf));
    if tempf<0
        tempf=tempf+pi;
    end
    phiit(j+1)=tempf+nte*pi+phiit(j);
end
%lamtau=phiit(i)/(a(i)*2^(i-1));
%un=
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
phin(:,i)=[(2.^(n-1)).*a(i)*h*ny+phiit(i);(2.^n)*a(i)*h*ny+phiit(i)];
while i > 1
    i = i - 1;
    phin(:,i) = phin(:,i+1);
    %if ~isempty(in)
      phin(:,i) = 0.5*(asin(c(i+1).*sin(rem(phin(:,i+1),2*pi))./a(i+1)) + phin(:,i+1));
    %end
end
sn = sin(rem(phin(:,1),2*pi));
cn = cos(rem(phin(:,1),2*pi));
dn = sqrt(1 - m .* sn.^2);

% special case m = 1 
m1 = find(m==1);
sn(m1) = tanh([h*ny/2+phiit(m1)/((2.^(n)).*a(i));h*ny+phiit(i)/((2.^(n)).*a(i))]);
cn(m1) = sech([h*ny/2+phiit(m1)/((2.^(n)).*a(i));h*ny+phiit(i)/((2.^(n)).*a(i))]);
dn(m1) = sech([h*ny/2+phiit(m1)/((2.^(n)).*a(i));h*ny+phiit(i)/((2.^(n)).*a(i))]);


