with(plots):

k1:=1/sqrt(2);

fsolve(JacobiSN(u,k1)=1,u);

plot({JacobiSN(u,k1),JacobiCN(u,k1),JacobiDN(u,k1)},
> u=-1.854074677..3*1.854074677);

delau:=[-a*(1-epsilon)*(EllipticK(kk)+
> EllipticF(JacobiSN(u,kk),kk))-a*(1+epsilon)*
> (EllipticE(kk)+EllipticE(JacobiSN(u,kk),kk)),a*(1+epsilon)*
> JacobiDN(u,kk)*cos(v),a*(1+epsilon)*JacobiDN(u,kk)*sin(v)];

with(plots):

plotund:=proc(a,eps,theta,phi)
> local kk,uu;
> kk:=2*sqrt(eps)/(1+eps);
> uu:=fsolve(JacobiSN(u,kk)=1,u);
> plot3d(subs({u=u(t),v=v(t)},[-a*(1-eps)*(EllipticK(kk)+
> EllipticF(JacobiSN(u,kk),kk))-a*(1+eps)*(EllipticE(kk)+
> EllipticE(JacobiSN(u,kk),kk)),a*(1+eps)*JacobiDN(u,kk)*
> cos(v),a*(1+eps)*JacobiDN(u,kk)*sin(v)]),u=-uu..uu,v=0..2*Pi,
> scaling=constrained,axes=boxed,orientation=[theta,phi],
> shading=XY,lightmodel=light3);
> end:

plotund(1,0.7,-79,70);

plotund(1,0.3,-65,70);

EFG(subs(kk=2*sqrt(epsilon)/(1+epsilon),delau));