function E=ellipticE(phi,k)
% E(phi,k)=integrate(sqrt(1-k sin(u)^2),u,[0,phi])
% Usage:: E=ellipticE(phi,k)

% revision history:
% 07/18/04 Mark D. Shattuck <mds> ellipjj.m
%          taken from Numerical Recipies 2nd Ed.
% 07/18/04 mds k->k^2 for NR notation

%% main

s=sin(phi);
c2=cos(phi).^2;
u=s.*sqrt(k);
v=(1-u).*(1+u);
E=s.*rf(c2,v,1)-1/3*k.*s.^3.*rd(c2,v,1);
