clear
close

load om_dat;
omega=om_dat(1);

load vzn_dat;
z=vzn_dat(1)
n=vzn_dat(2)+1
nl=vzn_dat(3);
for l=1:nl-1,
  thk(l)=vzn_dat(3+l)/1000;
end
load vr3_dat;
load vr4_dat;
nv=size(vr3_dat,1)
md=nv/n

for m=1:md,
  for j=1:n,
    h(j)=vr3_dat(n*(m-1)+j);
    v(j)=vr4_dat(n*(m-1)+j);
  end
  vmax=max(h);
  vmax2=max(v);
  if vmax2 > vmax
    vmax=vmax2;
  end
  vmin=min(h);
  vmin2=min(v);
  if vmin2 < vmin
    vmin=vmin2;
  end
  if -vmin > vmax
    vmax=-vmin;
  end
  for j=1:n,
    v1(j,m)=h(j)*0.45/vmax+m-1;
    v2(j,m)=v(j)*0.45/vmax+m-1;
  end
end

zd=z/(n-1);
for j=1:n,
  dz(j)=(-zd*(j-1))/1000;
end

plot(v1,dz,'-',v2,dz,'--')
hold on
x(1)=-0.5;
x(2)=md-0.5;
y(1)=0;
y(2)=0;
for l=1:nl-1,
  y(1)=y(1)-thk(l);
  y(2)=y(2)-thk(l);
  plot(x,y,':')
end
hold off
axis([-0.5 md-0.5 -z/1000 0])
xlabel('Mode (Solid: Horizontal, Dashed: Vertical)')
ylabel('Depth (km)')
ttl=['Stress Vectors of Rayleigh Wave at = ',num2str(omega/6.283185307), ' (Hz)'];
;
title(ttl)


print vrs.ps




