
function rayleigh_waves

clear all
close all

f = 0.5;
T = 1/f
omega = 2*pi*f
vel = 3000

ncol = 7

% Set up array of depths
x = [0:300:15000]
z = [0:-100:-2000]

Ax = 200;
Az = 200;

nx = length(x);
nz = length(z);
k = 2*pi*f/vel;

lamda = 750;

%Start time step loop
n=40;
M=moviein(n);
for it=1:n
  t=T*(it)/n;
  for ix=1:nx
    dx(ix) = 0.001*(x(ix)+Ax*exp(z(1)/lamda)*sin(k*x(ix)-omega*t));
    dz(ix) = 0.001*(z(1)+Az*exp(z(1)/lamda)*cos(k*x(ix)-omega*t));
  end
  
  plot(dx,dz,'o-')
  hold on
  plot(dx(20),dz(20),'or')
  axis([-0.2,15.2,-2.05,1])
  title('Rayleigh wave   f=0.5 Hz')
  ylabel ('depth(km)')
  xlabel ('(km)')
   
   
  for iz=2:ncol-1
    for ix=1:nx
      dx(ix) = 0.001*(x(ix)+Ax*exp(z(iz)/lamda)*sin(k*x(ix)-omega*t));
      dz(ix) = 0.001*(z(iz)+Az*exp(z(iz)/lamda)*cos(k*x(ix)-omega*t));
    end
    plot(dx,dz,'o')
  end
    
  for ix=1:nx
    dx(ix) = 0.001*(x(ix)+Ax*exp(z(ncol)/lamda)*sin(k*x(ix)-omega*t));
    dz(ix) = 0.001*(z(ncol)+Az*exp(z(ncol)/lamda)*cos(k*x(ix)-omega*t));
  end
  plot(dx,dz,'ro')
    
  for iz=ncol+1:nz
    for ix=1:nx
      dx(ix) = 0.001*(x(ix)+Ax*exp(z(iz)/lamda)*sin(k*x(ix)-omega*t));
      dz(ix) = 0.001*(z(iz)+Az*exp(z(iz)/lamda)*cos(k*x(ix)-omega*t));
    end
    plot(dx,dz,'o')
  end    
 
  hold off
  M(:,it)=getframe;
end
  
movie(M,10)
  
  
print -dpsc  rayleigh-wave.ps