clear
del=input('Enter delta');
wave=input('enter wave type (kel,ros,mrg,ing)','s');
k=input('enter wavenumber k');
n=input('Enter mode number n (max 3; not used for "kel")');
if wave == 'kel'
   n=-1;
end   
if wave == 'mrg'
   n=0;
end   
%
pi=acos(-1);
if wave == 'kel'
   w=k;
elseif wave == 'mrg'
   w=0.5.*(k-sqrt(k^2+4.*del));
elseif wave == 'ros'
   w=-k./(2.*n+1+k^2./del);
else
   w=sqrt(k.^2+del.*(2.*n+1));
end
%
y=[-2:0.2:2];
x=[0:0.1.*pi:pi];
ys=max(size(y));
xs=max(size(x));
ey=exp(-0.5.*y.^2);
h(1,1:ys)=1;
h(2,1:ys)=2.*y;
h(3,1:ys)=4.*y.^2-2;
h(4,1:ys)=8.*y.^3-12.*y;
%
hp(1,1:ys)=0;
hp(2,1:ys)=2;
hp(3,1:ys)=8.*y;
hp(4,1:ys)=24.*y.^2-12;
if wave == 'kel'
   for ix=1:xs,
      for iy=1:ys,
         u(ix,iy)=ey(iy).*cos(x(ix));
         s(ix,iy)=-u(ix,iy);
         wv(ix,iy)=k.*ey(iy).*sin(x(ix));
         wvp(ix,iy)=max(wv(ix,iy),0);
         wvn(ix,iy)=min(wv(ix,iy),0);
         v(ix,iy)=0;
      end
   end
else
   for ix=1:xs,
      for iy=1:ys,
         v(ix,iy)=ey(iy).*h(n+1,iy).*cos(k.*x(ix));
         u(ix,iy)=-ey(iy).*(k.*hp(n+1,iy)-w.*y(iy).*h(n+1,iy)).*sin(x(ix))./(k.^2-w.^2);
         s(ix,iy)=-w.*u(ix,iy)./k-y(iy).*h(n+1,iy).*ey(iy).*sin(x(ix))./k;
         wv(ix,iy)=ey(iy).*(w.*w.*hp(n+1,iy)-w.*y(iy).*k.*h(n+1,iy)).*cos(x(ix))./(k.^2-w.^2);
         wvp(ix,iy)=max(wv(ix,iy),0);
         wvn(ix,iy)=min(wv(ix,iy),0);
      end
   end
end
%
x=transpose(x);
y=transpose(y);
u=transpose(u);
v=transpose(v);
s=transpose(s);
wvp=transpose(wvp);
wvn=transpose(wvn);
j=pcolor(x,y,s);
shading interp
set(gca,'fontweight','bold')
hold on
g=quiver(x,y,u,v);
r=contour(x,y,wvp,'r');
q=contour(x,y,wvn,'b');
hold off
xlabel('kX','fontweight','bold')
ylabel('Y','fontweight','bold')
titstring=['s, u, v, and w, ' wave ',  n= ' num2str(n) '  k= ' num2str(k) '  delta= ' num2str(del)];
title(titstring,'fontweight','bold')