function cnoidal2(H,Period,d,z)
% 2nd Order Cnoidal Wave Theory
% Developed after ACES Technical Reference Material
close;
ep=H/d;
g=9.80;  rho=1020;
%	Complete elliptic integral of First and Second Kind
%Determine the value of k such that T=Period
kk=[0.0001:.01:.7 .701:.001:.999];
TT=elipp(kk,H,d,ep);
k=interp1(TT,kk,Period);
T=elipp(k,H,d,ep);
kp=sqrt(1-k.^2); lambda=kp.^2/k.^2; m=k.^2;
[K,E]=ellipke(m);
mu=E/(k.^2.*K);
%	Computation of Celerity
C0=1;
C1=(1+2*lambda-3*mu)/2;
C2=(-6-16*lambda+5*mu-16*lambda^2+10*lambda*mu+15*mu^2)/40;
c=sqrt(g*d)*(C0+ep*C1+ep^2*C2);
L=c*T;
x=0:.01:1; 
%	z=0:d/100:d;
t=0;	%	t=0:.1:1
%	Assume that all is at time 0 temporally 
theta=2*K*(x-t);
%	Computation of water surface elevation
A0=ep*(lambda-mu)+ep^2*((-2*lambda+mu-2*lambda^2+2*lambda*mu)/4);
A2=3/4*ep^2;
A1=ep-A2;
%	Compute Jacobian Elliptic Integral
[sn,cn,dn] = ELLIPJ(theta,m);
csd=sn.*cn.*dn;
eta=d*(A0+A1*cn.^2+A2*cn.^4);
%	Compute velocities
B20=-ep^2;
B00=ep*(lambda-mu)-B20*(lambda-mu-2*lambda^2+2*mu^2)/4;
B10=ep-B20*(1-6*lambda+2*mu)/4;
B01=-3*lambda*B20/2;
B11=-3*B20*(1-lambda);
B21=9*B20/2;
%z=input('Depth for desired velocities:')
if z<-d, display('z cannot be below mudline\n'),end
u=sqrt(g*d)*((B00+B10*cn.^2+B20*cn.^4)...
   -((z+d)/d).^2.*(B01+B11*cn.^2+B21*cn.^4)/2);
w=sqrt(g*d)*4*K*d*csd/L.*(((z+d)/d)*(B10+2*B20*cn.^2)...
   -((z+d)/d)^3*(B11+2*B21*cn.^2)/6);
%	Compute pressure
P0=3/2;
P1=(1+2*lambda-3*mu)/2;
P2=(-1-16*lambda+15*mu-16*lambda^2+30*lambda*mu)/40;
Pb=rho*g*d*(P0+ep*P1+ep^2*P2);
%	p=Pb-rho*((u-c)^2+w^2)/2-rho*g*(z+d);
%	Compute vertical velocity
for jj=1:size(u')
   if z>eta(jj)
      u(jj)=NaN;w(jj)=NaN;
   end
end
airy=H*cos(x*2*pi)/2;
subplot (2,1,1),plot(x,eta),title('Water Surface'),
			ylabel('Elevation(m)'),grid;hold on;
subplot (2,1,1), plot (x,airy,'r--');
subplot (2,1,2), plot (x,u),title('U and V'),
			ylabel('Velocity (m/s)'),grid; hold on;
subplot (2,1,2), plot (x,w,'g--'),xlabel('Horizontal Position');
%_____________________________________________
function T = elipp(k,H,d,ep)
%Determine T from k
m=k.^2;g=9.8;
[K,E]=ellipke(m);
mu=E/(k.^2.*K);
kp=sqrt(1-k.^2);
lambda=kp.^2/k.^2;
T=sqrt((16*k.^2.*K.^2*d^2)/(g*H*(1-ep.*(1+2*lambda)/4)));