% intwavesrho.m
% b klinger, 2/16/04

iwhich=input('timestep: 1 for short, 2 for long, 3 for very long: ');
xtot=30;
ztot=30;
N=2*pi/600;
if iwhich==1
   dt=30;
elseif iwhich==2
   dt=120;
elseif iwhich==3
   dt=480;
end
nt=30;
%nt=90;
k0=2*pi/5;
m0=.5*2*pi/5;
alpha=.1;
dk=alpha*k0;
delk=.1*k0;
dm=sqrt((1+1^2)-(1+alpha)^2)-1;
delm=(m0/k0)*delk;

x=[0:.01:1]*xtot;
z=[0:.01:1]*ztot;
nx=length(x);
nz=length(z);

k1=k0+dk+delk;
k2=k0+dk-delk;
k3=k0-dk+delk;
k4=k0-dk-delk;

m1=m0+dm+delm;
m2=m0+dm-delm;
m3=m0-dm+delm;
m4=m0-dm-delm;

omeg1=k1*N/sqrt(k1^2+m1^2);
omeg2=k2*N/sqrt(k2^2+m2^2);
omeg3=k3*N/sqrt(k3^2+m3^2);
omeg4=k4*N/sqrt(k4^2+m4^2);
per=2*pi./[omeg1 omeg2 omeg3 omeg4];
disp(per)

rho=zeros(nz,nx,nt);
w=zeros(nz,nx,nt);
for itime=1:nt
   t=dt*itime;
   for iz=1:nz
      sin1=sin(k1*x+m1*z(iz)-omeg1*t)/omeg1;
      sin2=sin(k2*x+m2*z(iz)-omeg2*t)/omeg2;
      sin3=sin(k3*x+m3*z(iz)-omeg3*t)/omeg3;
      sin4=sin(k4*x+m4*z(iz)-omeg4*t)/omeg4;
      rho(iz,:,itime)=(iz/nz)-(.02/(4*125))*(sin1+sin2+sin3+sin4);
      cos1=cos(k1*x+m1*z(iz)-omeg1*t);
      cos2=cos(k2*x+m2*z(iz)-omeg2*t);
      cos3=cos(k3*x+m3*z(iz)-omeg3*t);
      cos4=cos(k4*x+m4*z(iz)-omeg4*t);
      w(iz,:,itime)=cos1+cos2+cos3+cos4;
   end
end


for itime=1:nt
   hold off
   [cs,han]=contour(x,z,rho(:,:,itime),0:.1:1,'k');
   for i=1:length(han); set(han(i),'linewidth',1.5); end
   hold on
   contour2(x,z,w(:,:,itime),-3.5:3.5,'brww')
   axis equal; axis([0 30 0 30])
   grid
   t=dt*itime;
   title(['time = ' num2str(t) ' s'])
   pause(.1)
end
