%
%  Plots Love Wave solution, written by Jeff Gu
%
%  b1=input('enter shear velocity of layer1 [km/s]:');
%  r1=input('enter density of layer1 [kg/m_cube]:');
%% shear modulus, mu = b1*b1*10^6*r1*1000
 % b2=input('enter shear velocity of layer2 [km/s]:');
 % r2=input('enter density of layer2 [kg/m_cube]:');
 % h=input('enter thickness of layer over halfspace (km):');
 % n=input('enter total number of cx points:');
 % w=input('enter frequency omega in herz:');
set(0,'DefaultAxesColorOrder',[1 0 0;0 1 0;0 0 1],...
      'DefaultAxesLineStyleOrder','-|--|:');
% layer velocities v1 v2
  b1=3.9;
  b2=4.6;
% density in kilos of kg/m^3
  r1=2.8;
  r2=3.3;
  h=40;
  n=1000;
% mu values computed based on velocity and density
  mu1=b1*b1*1000000*r1*1000.;
  mu2=b2*b2*1000000*r2*1000.;
% frequency in Hz
  f=0.2;

  w=f*2.0*3.1415926;   
  dcx = (b2-b1)/n;
% range of phase velocity cx where b1<cx<b2
  cx=b1:dcx:b2;
  
  zeta=[];
  lterm=[];
  rterm=[];
  for i=1:n+1
     zeta(i)=h.*(cx(i)*cx(i)/(b1*b1)-1)^(0.5)/cx(i);
     lterm(i)=tan(w*zeta(i));
     rterm(i)=(mu2/mu1)*h*(1.-cx(i)*cx(i)/(b2*b2))^(0.5)/(cx(i)*zeta(i))
  end
  figure;
  subplot(3,1,1)
  plot(cx,zeta','.')
  xlabel('Phase velocity Cx');
  ylabel('Zeta');
  ylim([0 max(zeta)]);
  subplot(3,1,2)
  plot(zeta, lterm, '-', 'linewidth', 2);
  hold on
  plot(zeta, rterm, 'k-', 'linewidth', 2);
  xlabel('Zeta');
  ylabel('Love Wave Soln');
  ylim([-10 20]);
  subplot(3,1,3)
  plot(cx, lterm, 'linewidth', 2);
  hold on;
  plot(cx, rterm, 'k-', 'linewidth', 2);
  xlim([b1 b2]);
  ylim([-20 20]);
  xlabel('Cx');
  ylabel('Love Wave Soln');
  
