function [sn,cn,dn]=ellipticJ(u,E)
% ellipticJ <EllipticJ(u,E)>
% Usage:: [sn,cn,dn]=ellipticJ(u,E)
%
% Converted ellipj for energy E

% revision history:
% 07/18/04 Mark D. Shattuck <mds> ellipjj.m
% 07/20/04 mds fix E==1

%% Main

sn=zeros(max([size(u);size(E)]));
cn=sn;
dn=sn;
i1=find(E<1);
i2=find(E>1);
i3=find(E==1);

if (length(u)>1 && length(E)>1)
    [sn(i1) cn(i1) dn(i1)]=ellipj(u(i1),E(i1));
    [sn(i3) cn(i3) dn(i3)]=ellipj(u(i3),E(i3));
    [sn(i2) cn(i2) dn(i2)]=ellipj(u(i2).*sqrt(E(i2)+eps),1./(E(i2)+eps));
elseif (length(E)>1)
    [sn(i1) cn(i1) dn(i1)]=ellipj(u,E(i1));
    [sn(i3) cn(i3) dn(i3)]=ellipj(u,E(i3));
    [sn(i2) cn(i2) dn(i2)]=ellipj(u.*sqrt(E(i2)+eps),1./(E(i2)+eps));
else
    if(E<=1)
        [sn cn dn]=ellipj(u,E);
    else
        [sn cn dn]=ellipj(u*sqrt(E+eps),1/(E+eps));
    end
    
end    
