clear
close

load wn_dat;
size(wn_dat)
load qxr_dat;
load qxi_dat;
load phsl_dat;
nl=size(phsl_dat,1);
load phsr_dat;
nr=size(phsr_dat,1);
load in_dat;
om=in_dat(1);
ny=in_dat(2);
for i=1:ny
  vp(i)=in_dat(3+2*(i-1));
  vs(i)=in_dat(4+2*(i-1));
end

xmax=max(wn_dat);
y(1)=min(qxr_dat)*1.1;
y(2)=max(qxr_dat)*1.1;
subplot(2,1,1)
plot(wn_dat,qxr_dat,'-')
hold on
for i=1:nl
  x(1)=om/phsl_dat(i);
  x(2)=om/phsl_dat(i);
  plot(x,y,':')
end
for i=1:nr
  x(1)=om/phsr_dat(i);
  x(2)=om/phsr_dat(i);
  plot(x,y,'--')
end
for i=1:ny
  x(1)=om/vp(i);
  x(2)=om/vp(i);
  plot(x,y,'-.')
end
for i=1:ny
  x(1)=om/vs(i);
  x(2)=om/vs(i);
  plot(x,y,'-')
end
hold off
axis([0 xmax y(1) y(2)])
ylabel('Real Part of Qx')
title('Dotted: Love;  Dashed: Rayleigh;  Dashdot: Vp;  Solid: Vs')

y(1)=min(qxi_dat)*1.1;
y(2)=max(qxi_dat)*1.1;
if y(2) < -y(1)/1000
  y(2)=-y(1)*0.1
end
if -y(1) < y(2)/1000
  y(1)=-y(2)*0.1
end
subplot(2,1,2)
plot(wn_dat,qxi_dat,'-')
hold on
for i=1:nl
  x(1)=om/phsl_dat(i);
  x(2)=om/phsl_dat(i);
  plot(x,y,':')
end
for i=1:nr
  x(1)=om/phsr_dat(i);
  x(2)=om/phsr_dat(i);
  plot(x,y,'--')
end
for i=1:ny
  x(1)=om/vp(i);
  x(2)=om/vp(i);
  plot(x,y,'-.')
end
for i=1:ny
  x(1)=om/vs(i);
  x(2)=om/vs(i);
  plot(x,y,'-')
end
hold off
axis([0 xmax y(1) y(2)])
ylabel('Imaginary Part of Qx')
xlabel('Wavenumber')

print qx1.ps




