%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This script computes dispersion properties and energies for 
% pseudo-Rayleigh waves and other modes supported by a "Pekeris" 
% wave guide. This is a layer of water over a solid half-space.
% For Rayleigh modes this was first described by Stoneley (1926). 
% All modes are reviewed by Ardhuin and Herbers (JFM 2013)
%
% For storing the G transfer function: set storeGall = 1
% For plotting vertical mode structures: set doplot = 1
%                           and iplot1 < iplot2  (indices for which plots 
%                                                 will be made)
%
% Last modified: 10 Mar 2013
%
% History     :    Jun-2011 creation
%             : 10-Jul-2011 calculation of coupling coefficients 
%             : 19-Aug-2011 verification of energy integrals, add output file
%             : 30-Dec-2011 additional calculations for acoustic noise
%             : 31-Dec-2011 refinement for k estimation for mode 0
%             : 10-Mar-2013 adding transfer function for bottom pressure
%             : 20-Mar-2013 clean up 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script uses: 
%heading_sphere.m
%sismo_all_Pak135.txt
%sismo_Pekeris_response.m
%sismo_Rayleigh_disp_left.m
%sismo_Rayleigh_disp.m
%sismo_stoneley_all.m
%sismo_Stoneley_det_cof.m
%sismo_Stoneley_det.m
%sismo_Stoneley_disp2.m
%sismo_Stoneley_disp.m
%
% It also uses travel time data: sismo_all_Pak135.dat
%
% Defines the parameters of the 2 layers: values used by LH 1950
%
rhow =1000;                % density of water in kg/m^3
rhos =2.5*rhow;            % density of rock  in kg/m^3
betas=3000; %sqrt(10.E6);  %2800;     % shear wave velocity of rock in m/s
alphaw=1500; %sqrt(2.2E6);
%alphaw=1400; %sqrt(2.2E6);


%
% Alternative values ... used by Ardhuin et al. (JGR 2012)
%
alphaw=1470;
rhow =1026;                % density of water in kg/m^3
rhos =3300;                % density of rock  in kg/m^3
betas=3000; %sqrt(10.E6);  %2800;     % shear wave velocity of rock in m/s

alphas=sqrt(3)*betas;      % Compression velocity in case of Poisson solid

%
% Sediment layer: if needed, for comparison without sediment
% theory in Ardhuin et al. (submitted to JASA)
%
alphal=1500;
rhol=1700;
betal=25;
hl=100;    % sediment layer thickness 

g=0;
Q=g/(2*alphaw);
R=g/(2*betas);
S=g/(2*alphas);

h=3000;                    % water depth

tag ='_3000_c3';           % suffix for output file names

nfreq=300;   % we have used nfreq=1500 in Ardhuin and Herbers 2013           
freqs=linspace(1E-7,1.0,nfreq);     % wave frequencies 
nmodes=10;                          % number of vertical modes
ksolution=zeros(nfreq,nmodes)+NaN;  % array with wavenumbers 


%
% Compute critical frequencies (at which a new mode appears)
%
loms=-sqrt((betas^2-alphaw^2)/(alphas^2-betas^2))*alphas/alphaw*rhos/rhow;
ls=(atan(loms)+pi*linspace(0,nmodes-1,nmodes))/h;
sigmas=alphaw*ls/(sqrt(1-alphaw^2/betas^2));

ks=sigmas/betas;
Gprimks=ks.*0;

nks=50000;


storeGall=1 

if (storeGall == 1) 
  Gall2D=zeros(nfreq,nks);
  Gall2Dl=zeros(nfreq,nks);
end 

kk=linspace(1E-12,1E-2,nks);
dk=kk(2)-kk(1);

    
for ij=3:nmodes;
    sigma=sigmas(ij-1);
    k=ks(ij-1);
    fac=1+dk/k;
    [detpk,Gpk]=sismo_Stoneley_det(k*fac,sigma,rhow,rhos,alphaw,alphas,h,betas);
    Gprimk=(Gpk*sqrt((k*fac)^2-k^2));
    Gprimks(ij-1)=Gprimk;
end

% 10-Mar-2013: adding pressure transfer functions in 13 and 14
kdisp=zeros(nfreq,nmodes,15)+NaN;

PScoefs=zeros(nfreq,8)+NaN;
PSdir2=zeros(nfreq,91)+NaN;

iplot1=1;
iplot2=nfreq;
doplot=0;


    
for i=1:nfreq
       msg=sprintf('Computing response for frequency %f ,  %d out of  %d %d %d \n', ... 
      freqs(i),i,nfreq);
    fprintf('%s \n',msg)

    [i freqs(i)]
    sigma=2*2*pi*freqs(i);
    k01=sigma/(0.9*betas);
%
% finds wavenumber for the Rayleigh wave without water
%
    kR = fzero(@(k) sismo_Rayleigh_disp(k,sigma,alphas,betas),k01);
    ksolution(i,1)=kR;
%
% finds where the left hand side changes sign
%
    
    left=sismo_Rayleigh_disp_left(kk,sigma,alphaw,h);
    l=sqrt(sigma.^2./alphaw.^2-kk.^2);
    m=conj(sqrt(kk.^2-sigma.^2./alphas.^2));
    n=conj(sqrt(kk.^2-sigma.^2./betas.^2));
    bracket=(sigma.^2-2.*kk.^2.*betas.^2);
    right=(rhos/rhow).*(l./m).*(4.*betas.^4.*kk.^2.*(m.*n)-bracket.^2)./(sigma.^4);
    rightim=right;

    qp=rhow*sigma*m*sigma^2./bracket;
    rp=l./sigma.*((4.*betas.^4.*kk.^2.*(m.*n))./bracket-bracket);
    
    IR=find(imag(right)==0);
    
    changesign=ones(nks,1);
    changesign(2:nks)=real(left(2:nks)).*real(left(1:nks-1));
    I=find(changesign < 0 & abs(left') > 8 & imag(right')==0);
    nintervals=length(I)+1;
    
    intervals=zeros(nintervals,2);
% The first lower bound is the wavenumber for which n=0;    
    kmin=sigma/betas;
    intervals(1,1)=kmin;
    
    for ji=1:nintervals-1
        intervals(ji,2)=kk(I(ji)-1);
        intervals(ji+1,1)=kk(I(ji));
    end
    J=find(imag(left)>0) ;
% The last upper bound is the wavenumber for which l=0;    
    kmax=sigma/alphaw;
    intervals(nintervals,2)=kmax; %kk(J(1)-1);
%
%  gets the wavenumbers of the different modes: finds local minima in
%  function, starting from the last interval (this is mode 0)
%
    j0=1;
    for ji=0:min([nintervals-1 nmodes-1])
        if (sigma > sigmas(ji+1)) 
           k = fminbnd(@(k) sismo_Stoneley_disp2(k,sigma,rhow,rhos,alphaw,alphas,h,betas),... 
              intervals(nintervals-ji,1),intervals(nintervals-ji,2),optimset('TolX',1e-12,'Display','off'));
          

        
%
% Checks that this minimum corresponds to a true zero 
%
            if imag(sismo_Stoneley_disp(k,sigma,rhow,rhos,alphaw,alphas,h,betas)) == 0 & length(k) > 0
                ksolution(i,ji+2)=k;
                j0=j0+1;
            end
       end
    end
    
%
% prepares for a plot
%
    if (doplot == 1 & i >= iplot1 & i <= iplot2)
    I=find(abs(imag(right)) > abs(0.01*real(right)));
    right(I)=NaN;

    figure(30)
    clf
    plot(kk,real(left),'k+',kk,imag(left),'gx',kk(IR),right(IR),'ro')
    axis([5E-4 2E-3 -5 5])
    set(gca,'FontSize',16);
    xlabel('k (rad/m)')
    legend('Re(tan(lh))','Im(tan(lh))','RHS')
    grid
    hold on
    plot([kR kR],[-5,5],'g-','Linewidth',2)
    for ji=0:min([nintervals-1 nmodes-1])
    plot([ksolution(i,ji+2) ksolution(i,ji+2)],[-5,5],'b-','Linewidth',2)
    end
    hold off
    pause
    end
    
    IP=find(kk < sigma/alphaw);
    
    IP2=find(kk < sigma/alphas);
    IP3=find(kk < sigma/alphaw  & kk > sigma/betas);
    IP4=find(kk > sigma/alphas  & kk <= sigma/betas);
    
    [detall,Gall]=sismo_Stoneley_det(kk,sigma,rhow,rhos,alphaw,alphas,h,betas);
    betal=1000;
    
    % Remove comment to also compute result with sediment layer
    %Glay=sismo_2layer_det(kk(IP),sigma,rhow,rhos,alphaw,alphas,h,betas, ...
    %                               rhol,alphal,betal,hl);
                               
    
    [detall,Gall,GCall,GDall] = sismo_Stoneley_det_cof(kk,sigma,rhow,rhos,alphaw,alphas,h,betas);
    if (storeGall == 1) 
       %[detall,Gall,GCall,GDall] = sismo_Stoneley_det_cof(kk(IP),sigma,rhow,rhos,alphaw,alphas,h,betas);
       l=sqrt(sigma.^2./alphaw.^2-kk.^2);

       % Transfer function with sediment layer
       %[Glay,GClay,GDlay]=sismo_2layer_det(kk,sigma,rhow,rhos,alphaw,alphas,h,betas, ...
       %                            rhol,alphal,betal,hl);
       Gall2D(i,:)=Gall;  
       
       
       % transfer to bottom pressure
       % l=sqrt(sigma.^2./alphaw.^2-kk.^2);   vertical wavenumber
       %Gall2D(i,:)=(GCall+GDall)*rhow*sigma; %all.*(rhow.*sigma.^2.*sin(l.*z)./(l.*cos(l.*h)));  % transfer to bottom pressure
       %Gall2Dl(i,:)=(GClay+GDlay)*rhow*sigma;   % transfer to bottom pressure
    end
    Gallnosing=Gall;
    dir=atan2(kk(IP2),abs(m(IP2)))*180./pi;
    if (length(IP2) > 2)
    dir2=linspace(0,90,91);

    dthet=dir(3:end)-dir(1:end-2);
    dthet=[dthet(1) dthet dthet(end)];
    
    % division by sin (phi) to have something per unit solid angle
    PSdir=rhos^2*betas^4.*(2.*abs(l(IP2).*m(IP2))./(abs(detall(IP2)))).^2./(alphas.*sigma).*dk./dthet.*180./pi;
    PSdir2(i,:)=interp1(dir,PSdir,dir2);
    end
    PScoefs(i,1)=2*pi*rhos^2*betas^4*sum(2.*abs(l(IP2).*m(IP2)./(sigma.*abs(detall(IP2)))).^2.*kk(IP2))*dk;
    
    PScoefs(i,1)=2*pi*rhos^2*betas^4*sum(2.*abs(l(IP2).*m(IP2)./(sigma.*abs(detall(IP2)))).^2.*kk(IP2))*dk;
    PScoefs(i,7)=2*pi*rhos^2*betas^4*sum(2.*abs(l(IP2).*m(IP2).*kk(IP2).^2  ...
                   ./((n(IP2).^2+kk(IP2).^2).*sigma.*abs(detall(IP2)))).^2.*kk(IP2))*dk;
    PScoefs(i,2)=2*pi*rhos^2*betas^4*sum(abs(Gall(IP2)).^2.*kk(IP2))*dk;

    PScoefs(i,5)=2*pi*rhow^2*alphaw^4*sum(abs(GCall(IP2)+GDall(IP2)).^2.*kk(IP2))*dk;
    
%
% Energy of Rayleigh waves
%
    kdisp(i,1,1)=(4.*pi.*freqs(i))/kR;
    ij=1;
    kdisp(i,ij,6)=kR;
    k=kR;
    m=sqrt(k.^2-sigma.^2./alphas.^2);
    n=sqrt(k.^2-sigma.^2./betas.^2);

    DeltatoA=m*sigma/(2*k^2*betas^2-sigma^2);
    BtoA=2*k*m/(n^2+k^2);
 
    Ecs=0.5*rhos*0.5*((k.*BtoA+m).^2/(2*m)+(n.*BtoA+k).^2/(2*n))/(DeltatoA^2);
%
%  without water the energy is only in the crust: 
%
    kdisp(i,1,2)=Ecs;
% 
%  Group velocity: Rayleigh waves (without the water layer) are not dispersive 
%
    kdisp(i,1,7)=kdisp(i,1,1);
    
  
    for ij=2:nmodes
      k=ksolution(i,ij);
      if (ij==2) 
%          if(2*freqs(i).*h./alphaw > 1.68) 
%              k=NaN;
%          end
      end
      if (isfinite(k)==1) 
 
         l=sqrt(sigma.^2./alphaw.^2-k.^2);
         m=sqrt(k.^2-sigma.^2./alphas.^2);
         n=sqrt(k.^2-sigma.^2./betas.^2);

         B2toA=m*sigma^2/((2*k^2*betas^2-sigma^2)*2*l*cos(l*h));
         DeltatoA=m*sigma/(2*k^2*betas^2-sigma^2);
         B2toDelta=sigma/(2*l*cos(l*h));
         BtoA=2*k*m/(n^2+k^2);
 
%
% The kinetic energy in the crust is the mean sum of the vertical velocity 
% variance w^2 = (d phi / dz + d psi / dx)^2   + u^2= (d phi / dx - d psi / dz)^2 
%
% variance w^2=(d phi / dz)^2   +    u^2=(d phi / dx)^2 
%
         Ecs=0.5*rhos*0.5*((k.*BtoA+m).^2/(2*m)+(n.*BtoA+k).^2/(2*n))/(DeltatoA^2);



%
% The kinetic energy in the water is rhow times the mean of half the velocity 
% variance w^2=(d phi / dz)^2   +    u^2=(d phi / dx)^2 
%
%      var(w)=0.5*[delta*2*l*cos(lz)*B2]^2   
%      var(u)=0.5*[delta*2*k*sin(lz)*B2]^2

         Ecww=0.5*rhow*0.5*(2*l)^2*(h/2+(sin(2*l*h)/(l*4)));
         Ecwu=0.5*rhow*0.5*(2*k)^2*(h/2-(sin(2*l*h)/(l*4)));
         Ecw=(Ecwu+Ecww).*(B2toDelta).^2;
        
       % [k l m n Ecs Ecww Ecwu Ecw B2toDelta]
                        
      if (doplot == 1 & i >= iplot1 & i <= iplot2)
         figure(40+ij)
         clf
         zz=linspace(-h*3,-h,600);
         w=( (m.*exp((zz+h).*m)) ...
                        -(k.*exp((zz+h).*n).*BtoA))./DeltatoA;   
         u=( (k.*exp((zz+h).*m)) ...
                        +(n.*exp((zz+h).*n).*BtoA))./DeltatoA;   
         plot(w,zz,'r--',u,zz,'k--','LineWidth',2);
         hold on;

         zz=linspace(-h,0,200);
         w=-2*l*cos(l.*zz).*B2toDelta;   
         u=2*k*sin(l.*zz).*B2toDelta;   
         plot(w,zz,'r-',u,zz,'k-','LineWidth',2);
         hold off;
         pause
       end

        kdisp(i,ij,1)=(4.*pi.*freqs(i))/k;
        kdisp(i,ij,2)=Ecw+Ecs;
        kdisp(i,ij,3)=l;
        kdisp(i,ij,4)=m;
        kdisp(i,ij,5)=n;
        kdisp(i,ij,6)=k;
        
%
% Group velocity =  d omega / d k
%
        if i > 2 
            kdisp(i-1,ij,7)=4*pi*(freqs(i)-freqs(i-2))./(kdisp(i,ij,6)-kdisp(i-2,ij,6));
        end
        kdisp(i,ij,11)=kdisp(i,ij,2).*kdisp(i,ij,7)./kdisp(i,ij,6)./(kdisp(i,1,2).*kdisp(i,1,7)./kdisp(i,1,6));
        kdisp(i,ij,12)=sigma/(2*l*cos(l*h));


      end
    end
%end


%
% Now computes G' for growth rates
%
    sigma=2*2*pi*freqs(i);
    kdisp(i,1,8)=0.;
    kdisp(i,1,9)=0.;
    kdisp(i,1,14)=0.;
    sing0=0;
    for ij=2:nmodes
      fac=1.00005;
      k=kdisp(i,ij,6);
      if (isfinite(k)==1) 
           if (isfinite(kdisp(i,min([ij+1,nmodes]),6))==0) 
               fac=min([sqrt(k/kmin),fac]);
           end
           if (ij == 2 & (kmax-k)< k*(fac-1) ) 
               %fac=min([(kmax/k),fac]);
               [detp,Gm,GPm]=sismo_Stoneley_det(k,sigma*fac,rhow,rhos,alphaw,alphas,h,betas);

               [detmk,Gmk,GPmk]=sismo_Stoneley_det(k/fac,sigma,rhow,rhos,alphaw,alphas,h,betas);
               % for vertical ground displacement
               Gprimk=Gmk*((k/fac)^2-k^2);
               Gprime=Gm*((sigma*fac)^2-sigma^2);
               % now for pressure
               GPprimk=GPmk*((k/fac)^2-k^2);
               GPprime=GPm*((sigma*fac)^2-sigma^2);
               detpk=-detmk;
           else 
              [detp,Gp, GPp]=sismo_Stoneley_det(k,sigma*fac,rhow,rhos,alphaw,alphas,h,betas);
              [detm,Gm, GPm]=sismo_Stoneley_det(k,sigma/fac,rhow,rhos,alphaw,alphas,h,betas);

              [detpk,Gpk, GPpk]=sismo_Stoneley_det(k*fac,sigma,rhow,rhos,alphaw,alphas,h,betas);
              [detmk,Gmk, GPmk]=sismo_Stoneley_det(k/fac,sigma,rhow,rhos,alphaw,alphas,h,betas);
          
%
% Computes dG/dsigma  and dG/dk with finite differences
%
              Gprime=0.5*(Gp*((sigma*fac)^2-sigma^2)+Gm*((sigma/fac)^2-sigma^2));
              Gprimk=0.5*(Gpk*((k*fac)^2-k^2)+Gmk*((k/fac)^2-k^2));

              GPprime=0.5*(GPp*((sigma*fac)^2-sigma^2)+GPm*((sigma/fac)^2-sigma^2));
              GPprimk=0.5*(GPpk*((k*fac)^2-k^2)+GPmk*((k/fac)^2-k^2));
           end
%               
% case when k corresponds to kmin: singularity like 1/sqrt(k'-k)
%
           if (abs(k-kmin) < 0.3*dk & ij > 2) 
                    fac=1+dk/k;
                  [detpk,Gpk, GPpk]=sismo_Stoneley_det(k*fac,sigma,rhow,rhos,alphaw,alphas,h,betas);
                  [detmk,Gmk, GPmk]=sismo_Stoneley_det(k/fac,sigma,rhow,rhos,alphaw,alphas,h,betas);
                  Gprime=0.; 
                  Gprimk=(Gpk*sqrt((k*fac)^2-k^2));
                  GPprime=0.; 
                  GPprimk=(GPpk*sqrt((k*fac)^2-k^2));
        % Removes right-side of singularity
                  %Gallnosing(IP3)= Gallnosing(IP3)-Gprimk./sqrt(kk(IP3).^2-k^2);
        % Defines integral of right side 
                  sing0=Gprimk/sqrt(2*k);
        % Removes left-side of singularity
                  %Gprimk=(Gmk*sqrt((k/fac)^2-k^2));
   
                  %TEST        
    %              Gallnosing(IP4)= Gallnosing(IP4)-Gprimk./sqrt(kk(IP4).^2-k^2);
                  Gallnosing= Gallnosing-Gprimks(ij-1)./sqrt(((kk.^2-sigma.^2/betas.^2)));
                  IZE=find(abs(kk(IP)-k) < dk);
                  if (IP(IZE) > 1 & IP(max(IZE)) < max(IP) ) 
                      Gallnosing(IP(IZE))=Gallnosing(IP(max(IZE))+1);
                  end
                  
                else
%
% If the determinant changes sign: singularity is like 1/(k'-k)
%
                  if (imag(detpk)==0 & imag(detmk)==0)
                    if (real(detpk)*real(detmk) < 0) 
                      Gallnosing(IP3)= Gallnosing(IP3)-Gprimk./(kk(IP3).^2-k^2);
                      IZE=find(abs(kk(IP)-k) < 0.45*dk);
                      Gallnosing(IP(IZE))=Gallnosing(IP(IZE)-1);
                    else
                   
%               
% new estimation of k 
%
                     facnew=1+0.5*dk/k;
                     knew= fminbnd(@(k) sismo_Stoneley_disp2(k,sigma,rhow,rhos,alphaw,alphas,h,betas),... 
                     k/facnew,k*facnew,optimset('TolX',1e-12,'Display','off'));
                     [detpk,Gpk,GPpk]=sismo_Stoneley_det(knew*fac,sigma,rhow,rhos,alphaw,alphas,h,betas);
                     [detmk,Gmk,GPmk]=sismo_Stoneley_det(knew/fac,sigma,rhow,rhos,alphaw,alphas,h,betas);
                     if (imag(detpk)==0 & imag(detmk)==0)
                       if (real(detpk)*real(detmk) < 0) 
                         k=knew;
                         [detp,Gp,GPp]=sismo_Stoneley_det(k,sigma*fac,rhow,rhos,alphaw,alphas,h,betas);
                         [detm,Gm,GPm]=sismo_Stoneley_det(k,sigma/fac,rhow,rhos,alphaw,alphas,h,betas);
                         Gprime=0.5*(Gp*((sigma*fac)^2-sigma^2)+Gm*((sigma/fac)^2-sigma^2));
                         GPprime=0.5*(GPp*((sigma*fac)^2-sigma^2)+GPm*((sigma/fac)^2-sigma^2));
                         Gprimk=0.5*(Gpk*((k*fac)^2-k^2)+Gmk*((k/fac)^2-k^2));
                         Gallnosing= Gallnosing-Gprimk./(kk(IP).^2-k^2);
                         IZE=find(abs(kk(IP)-k) < 0.45*dk);
                         Gallnosing(IP(IZE))=Gallnosing(IP(IZE)-1);
                       end
                     end
                   end
                
              end
            end  
               
%  This is the T_delta(omega) of Hasselmann 1963
        kdisp(i,ij,10)=pi*Gprime^2/(2*sigma^2);
        kdisp(i,ij,8)=pi*Gprime^2/(2*sigma^2)*k;
%  This is an analogue tranfer function for bottom pressure 
        kdisp(i,ij,13)=pi*GPprime^2/(2*sigma^2)*k;
        

      else  %%% removes sigularity in S wave domain
         if (abs(sigma-sigmas(ij-1))/betas < 15*dk & ij > 2) 
            %              Gallnosing(IP4)= Gallnosing(IP4)-Gprimk./sqrt(kk(IP4).^2-k^2);
                  Gallnosing= Gallnosing-Gprimks(ij-2)./sqrt(((kk.^2-sigma.^2/betas.^2)));
                  IZE=find(abs(kk(IP)-k) < dk);
                  if (IP(IZE) > 1 & IP(max(IZE)) < max(IP) ) 
                      Gallnosing(IP(IZE))=Gallnosing(IP(max(IZE))+1);
                  end

         end
      end
    end
          
      if (doplot == 1 & i >= iplot1 & i <= iplot2)
         figure(34)
         clf
        plot(kk,abs(Gall),'k',kk,abs(Gallnosing),'r')
          axis([0  sigma/alphas*1.5 -4E-3 4.E-3])
          set(gca,'FontSize',16);
         xlabel('k (rad/m)')
         legend('|G|','G no sing')
         grid
         pause
       end
       PScoefs(i,3)=2*pi*rhos^2*betas^4*sum(abs(Gallnosing(IP)).^2.*kk(IP))*dk  ...
           + sing0.^2*pi/4*kmin;
       PScoefs(i,4)=2*pi*rhos^2*betas^4*sum(abs(Gallnosing(IP)).^2.*kk(IP))*dk  ...
           + sing0.^2*pi/4*kmin;
    l=sqrt(sigma.^2./alphaw.^2-kk.^2);
    m=conj(sqrt(kk.^2-sigma.^2./alphas.^2));
    n=conj(sqrt(kk.^2-sigma.^2./betas.^2));

        PScoefs(i,8)=PScoefs(i,7)+2*pi*rhos^2*betas^4*sum(2.*abs(l(IP4).*m(IP4).*kk(IP4).^2  ...
                   ./((n(IP4).^2+kk(IP4).^2).*sigma.*abs(detall(IP4)))).^2.*kk(IP4))*dk;
    
end

%
%  Finishes job for growth rates using group velocities ... 
%

for i=1:nfreq
    PScoefs(i,6)=0.;
    sigma=2*2*pi*freqs(i)';
    for ij=2:nmodes
      k=kdisp(i,ij,6);
      if (isfinite(k)==1) 
           kdisp(i,ij,8)=kdisp(i,ij,8)./(kdisp(i,ij,7).^2);
%  This is c_n 
           kdisp(i,ij,9)=sqrt((kdisp(i,ij,8)/sigma)*betas^5*rhos^2/(2*pi));
           
               l=sqrt(sigma.^2./alphaw.^2-k.^2);
               m=sqrt(k.^2-sigma.^2./alphas.^2);
               n=sqrt(k.^2-sigma.^2./betas.^2);
               z=-h;
               ilzph=sqrt(-1).*l.*(z+h);
               fac1=(sigma^2-2*k.^2.*betas^2);
               fac2=(4.*betas.^4.*k.^2.*(m.*n)./fac1-fac1);
               rp=l.*rhos.*fac2./sigma;
               qp=rhow.*m.*sigma.^3./fac1; 
               GC=(sqrt(-1).*rp-qp);
               GD=(sqrt(-1).*rp+qp);

               %delta2top2=abs(GC.*exp(ilzph)+GD.*exp(-1.*ilzph))^2*abs(k^2+n^2)^2/(abs(2*l*m*(k^2-n^2))^2);
               %PScoefs(i,6)=PScoefs(i,6)+kdisp(i,ij,9).^2.*delta2top2;

               DeltatoA=m*sigma/(2*k^2*betas^2-sigma^2);
               delta2top2=abs(rhow*sigma^2*sin(l*z)/(l*cos(l*h)))^2;
               PScoefs(i,6)=PScoefs(i,6)+kdisp(i,ij,9).^2.*delta2top2;
          
           % without the sin(lz) factor
           delta2top2=abs(rhow*sigma^2/(l*cos(l*h)))^2;
           kdisp(i,ij,15)=sqrt(delta2top2);

               
           %kdisp(i,ij,13)=kdisp(i,ij,13)./(kdisp(i,ij,7).^2);
           kdisp(i,ij,13)=kdisp(i,ij,8).*kdisp(i,ij,15);
           kdisp(i,ij,15)=sqrt(delta2top2)*sin(l*z)./(sigma.*betas*rhow);
            
           %  This is a similar c_n for bottom pressure
           kdisp(i,ij,14)=sqrt((kdisp(i,ij,13)/sigma^3)*betas^3*rhos^2/(2*pi.*rhow^2));
           if (isfinite(kdisp(i,ij,9))==1) 
               kdisp(i,1,9)=kdisp(i,1,9)+kdisp(i,ij,9).^2;
               kdisp(i,1,14)=kdisp(i,1,14)+kdisp(i,ij,14).^2;
           end

%
%  SDF(kx,ky)=  c_n^2 * omega * 2*pi /(beta^5 * rhos^2)     
%


      end
    end
    kdisp(i,1,9)=sqrt(kdisp(i,1,9));
    kdisp(i,1,14)=sqrt(kdisp(i,1,14));
end



%
% Save growth rates to ASCII file for all modes
%
fid2=fopen(['Rayleigh_source' tag '.txt'],'w');
for i=1:nfreq
    fndim=2*freqs(i).*h./alphaw;
    sigma=2*2*pi*freqs(i);
    ci=[];
    for ij=2:nmodes
      ci=[ci sqrt((kdisp(i,ij,8)/sigma)*betas^5*rhos^2/(2*pi))];
    end
    A=[fndim,ci ];
    fprintf(fid2,' %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f \n',transpose(A));
end
fclose(fid2);

fid2=fopen(['Rayleigh_source_PB' tag '.txt'],'w');
for i=1:nfreq
    fndim=2*freqs(i).*h./alphaw;
    sigma=2*2*pi*freqs(i);
    ci=[];
    for ij=2:nmodes
      ci=[ci sqrt((kdisp(i,ij,13)/sigma)*betas^5*rhos^2/(2*pi))];
    end
    A=[fndim,ci ];
    fprintf(fid2,' %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f \n',transpose(A));
end
fclose(fid2);

fid2=fopen(['Rayleigh_modes_k_' tag '.txt'],'w');
for i=1:nfreq
    fndim=2*freqs(i).*h./alphaw;
    sigma=2*2*pi*freqs(i);
    ci=[];
    for ij=2:nmodes
      ci=[ci kdisp(i,ij,6)];
    end
    A=[fndim,ci ];
    fprintf(fid2,' %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f %8.6f  \n',transpose(A));
end
fclose(fid2);
%
% Save results to ASCII file for mode 0
%
fid2=fopen(['Rayleigh1_speed_lin' tag '.txt'],'w');
for i=1:nfreq
    fndim=2*freqs(i).*h./alphaw;
    %ratio=kdisp(i,2,2).*kdisp(i,2,7)./(kdisp(i,1,2).*kdisp(i,1,7));
    ratio=kdisp(i,2,2).*kdisp(i,2,7)/(kdisp(i,1,2).*kdisp(i,1,7));
    if isfinite(ratio)==0.
        ratio=1;
    end
    A=[fndim,kdisp(i,2,1),kdisp(i,2,7),kdisp(i,1,9),kdisp(i,2,9),kdisp(i,2,10),ratio];
    %if (fndim < 1.37)      
            fprintf(fid2,' %8.6f %8.2f %8.2f %8.4f %8.4f %16.10e %8.6f \n',transpose(A));
    %end
end
fclose(fid2);
fid2=fopen(['Rayleigh2_speed_lin' tag '.txt'],'w');
I=find(isfinite(kdisp(:,3,1))==1);
J=I(2);
CgE_conv=(kdisp(J,2,2).*kdisp(J,2,7))*kdisp(J,2,2).*kdisp(J,2,7)/(kdisp(J,1,2).*kdisp(J,1,7));
ratio=zeros(nfreq,1);
for i=1:nfreq
    fndim=2*freqs(i).*h./alphaw;
    %ratio=kdisp(i,2,2).*kdisp(i,2,7)./(kdisp(i,1,2).*kdisp(i,1,7));
    ratio(i)=kdisp(i,3,2).*kdisp(i,3,7)/CgE_conv;
    if isfinite(ratio(i))==0.
        ratio(i)=1;
    end
    A=[fndim,kdisp(i,3,1),kdisp(i,3,7),kdisp(i,3,8),kdisp(i,3,9),kdisp(i,3,10),ratio(i)]
    A0=[fndim,kdisp(i,2,1),kdisp(i,2,7),0,0,0,ratio(i)];
    if (kdisp(i,3,9) > 0.04 | fndim > 1.5  )      
            fprintf(fid2,' %8.6f %8.2f %8.2f %8.4f %8.4f %16.10e %8.6f \n',transpose(A));
    else
            fprintf(fid2,' %8.6f %8.2f %8.2f %8.4f %8.4f %16.10e %8.6f \n',transpose(A0));
    end
        
end
fclose(fid2);
   
%  Plots Longuet-Higgins' coefficients c_n
figure(1)
clf
plot(2*freqs(:).*h./alphaw,abs(kdisp(:,2:nmodes,9)),'LineWidth',2)
set(gca,'FontSize',16);
xlabel('f_s h / \alpha_w')
ylabel('c_n^2')
hold on
plot(2*freqs(:).*h./alphaw,kdisp(:,1,9),'k--','LineWidth',2)
set(gca,'FontSize',16);
xlabel('f_s h / \alpha_w')
ylabel('coupling coefficients')
legend('mode 0','mode 1','mode 2','mode 3','mode 4','mode 5','mode 6','all')
grid on
fs0=1/6;
h0=100;
fs0*h/alphaw;
hold off


%  Plots the equivalent of Longuet-Higgins' coefficients c_n but for bottom
%  pressure
figure(2)
clf
plot(2*freqs(:).*h./alphaw,abs(kdisp(:,2:nmodes,9).*kdisp(:,2:nmodes,15)),'LineWidth',2)
set(gca,'FontSize',16);
hold on
plot(2*freqs(:).*h./alphaw,kdisp(:,1,14),'k--','LineWidth',2)
set(gca,'FontSize',16);
xlabel('f_s h / \alpha_w')
ylabel('coupling coefficients for p_b')
legend('mode 0','mode 1','mode 2','mode 3','mode 4','mode 5','mode 6','all')
grid on
fs0=1/6;
h0=100;
fs0*h/alphaw;
hold off



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%
if (storeGall == 1) 
nfp=floor(nfreq*400/1500);
%figure(3)
%clf
%pcolor(kk(1:5000).*1000,2*freqs(1:nfp),10.*log10(abs(Gall2Dl(1:nfp,1:5000))));shading flat;
%caxis([-40 40]);
%colorbar;
%set(gca,'FontSize',16);

figure(4)
pcolor(kk(1:5000).*1000,2*freqs(1:nfp),10.*log10(abs(Gall2D(1:nfp,1:5000))));shading flat;
caxis([-120 -40]);
colorbar;
set(gca,'FontSize',16);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(5)
  plot(2*freqs(:).*h./alphaw,kdisp(:,1:nmodes,1),'LineWidth',2)
set(gca,'FontSize',16);
xlabel('f_s h / \alpha_w')
ylabel('phase speed (m/s)')
%legend('R','mode 0','mode 1','mode 2','mode 3','mode 4')


figure(6)
  plot(2*freqs(:).*h./alphaw,kdisp(:,1:nmodes,7),'LineWidth',2)
set(gca,'FontSize',16);
xlabel('f_s h / \alpha_w')
ylabel('Group speed')
legend('R','mode 0','mode 1','mode 2','mode 3','mode 4','mode 5','mode 6')

figure(7)
  plot(2*freqs(:).*h./alphaw,atan(kdisp(:,2:nmodes,3)./kdisp(:,2:nmodes,6)).*180./pi,'LineWidth',2)
set(gca,'FontSize',16);
xlabel('f_s h / \alpha_w')
ylabel('sound propagation angle in water')
legend('mode 0','mode 1','mode 2','mode 3','mode 4')

figure(8)
plot(2*freqs(:).*h./alphaw,kdisp(:,2:nmodes,2)./repmat(kdisp(:,1,2),1,nmodes-1),'LineWidth',2)
set(gca,'FontSize',16);
xlabel('f_s h / \alpha_w')
ylabel('Energy with water / energy without')
legend('mode 0','mode 1','mode 2','mode 3','mode 4','mode 5','mode 6')
axis([0 3.6 0 8])

figure(9)
% expected amplification factor when going from sea to land 
plot(2*freqs(:).*h./alphaw,kdisp(:,2:nmodes,2).*kdisp(:,2:nmodes,7)./repmat(kdisp(:,1,2).*kdisp(:,1,7),1,nmodes-1),'LineWidth',2)
%plot(2*freqs(:).*h./alphaw,kdisp(:,2:nmodes,2)./repmat(kdisp(:,1,2),1,nmodes-1),'LineWidth',2)
set(gca,'FontSize',16);
xlabel('f_s h / \alpha_w')
ylabel('Shoaling coefficient')
legend('mode 0','mode 1','mode 2','mode 3','mode 4','mode 5','mode 6')
axis([0 3.6 0 20])

figure(10)
plot(2*freqs(:).*h./alphaw,kdisp(:,2:nmodes,3).*h./(pi),'LineWidth',2)
set(gca,'FontSize',16);
xlabel('f_s h / \alpha_w')
ylabel('lh / \pi')
legend('mode 0','mode 1','mode 2','mode 3','mode 4','mode 5','mode 6')
axis([0 3.6 0 8])

figure(11)
clf
set(gca,'FontSize',14)
%plot(2*freqs(2:nfreq).*h./alphaw,sqrt(PScoefs(2:nfreq,3)),'LineWidth',2)
hold on
%plot(2*freqs(2:nfreq).*h./alphaw,sqrt(PScoefs(2:nfreq,2)),'k-','LineWidth',2)
plot(2*freqs(2:nfreq).*h./alphaw,sqrt(PScoefs(2:nfreq,1)),'r-','LineWidth',2)
plot(2*freqs(2:nfreq).*h./alphaw,sqrt(PScoefs(2:nfreq,7)),'b','LineWidth',2)
plot(2*freqs(2:nfreq).*h./alphaw,sqrt(PScoefs(2:nfreq,8)),'k','LineWidth',2)
%legend('bound R + P +S','P + S','S','S all','P')
ylabel('coupling coefficients c_P or c_S')
legend('P','S for k < \alpha_s','S for all k')
xlabel('f_s h / \alpha_w')
hold off

fndim=2*freqs(:).*h./alphaw;
fid2=fopen(['P_source' tag '.txt'],'w');
for i=1:nfreq
    
    sigma=2*2*pi*freqs(i);
    A=[fndim(i) sqrt(PScoefs(i,1))  sqrt(PScoefs(i,7)) sqrt(PScoefs(i,8))];
    fprintf(fid2,' %8.6f %8.6f %8.6f %8.6f  \n',transpose(A));
end
fclose(fid2);

PSdir2(:,1)= PSdir2(:,2);
for i=2:nfreq
    I=find(isfinite(PSdir2(i,:))==0);
    if (length(I) > 0) 
        PSdir2(i,I)=PSdir2(i,I(1)-1);
    end
end

varphi=dir2.*pi./180;
mu2=(alphas./betas)^2;

%  Amplification factor for body waves as a function of frequency and
%  take-off angle
figure(12)
pcolor(fndim,dir2,10.*log10(PSdir2'));shading flat;colorbar;
set(gca,'LineWidth',2,'FontSize',16,'Position',[0.15 0.165312 0.931715 0.756978]);
xlabel('f_s h / \alpha_w')
ylabel('take-off angle (degrees)')
axis([0 3 0 90])
colorbar
caxis([-40 20])


%%%%%%%%%%%%%%%%%%%%%%%%%%% r and q %%%%%%%%%%%%%

figure(13)
%        kdisp(i,ij,3)=l;
%        kdisp(i,ij,4)=m;
%        kdisp(i,ij,5)=n;
%        kdisp(i,ij,6)=k;

qq=zeros(nfreq,nmodes);
rr=qq;
om2=(2*pi*freqs).^2';
for i=1:nmodes
    k2=kdisp(:,i,6).^2;
    qq(:,i)=rhow.*(2*pi*freqs').*kdisp(:,i,4).*om2./(om2-2.*k2.*betas.^2);
    rr(:,i)=rhos.*kdisp(:,i,3)./(2*pi*freqs').*(om2-2.*k2 ...
         -4.*betas.^4.*k2.*kdisp(:,i,4).*kdisp(:,i,5)./(om2-2.*k2.*betas^2));
end

plot(2*freqs(:).*h./alphaw,abs(rr./qq),'LineWidth',2)
set(gca,'FontSize',16);
xlabel('f_s h / \alpha_w')
ylabel('Effective impedance ratio r/q')
legend('mode 0','mode 1','mode 2','mode 3','mode 4')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Body wave propagation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% loads results from iaspei-tau : travel times for body waves 

Pdata=load('sismo_all_Pak135.txt');  % type 0 is Pg, 1 is Pn, 2 is P, 3 is PP. 
Pang=Pdata(:,5);               % take-off angle
Ptyp=Pdata(:,3);
Ptau=Pdata(:,4);               % travel time
Pdel=Pdata(:,1);               % angle at the center of the Earth
nP=length(Ptyp);
for i=2:nP
    if (Pdel(i)==Pdel(i-1) & Ptyp(i)==Ptyp(i-1))
        Ptyp(i)=Ptyp(i)+1;
    end
end
IP=find(Ptyp==3);

figure(14)
plot(Pdel(IP(1:end-1)),(Pang(IP(1:end-1))-Pang(IP(2:end)))./(Pdel(IP(1:end-1))-Pdel(IP(2:end))))

figure(15)
plot(Pdel(IP),Pang(IP))
i=351;


ddel=50/40000*360;



dOmega=zeros(length(IP),1);
for i=2:length(IP)-1
    if abs(Pang(IP(i+1))-Pang(IP(i-1))) < 2
    phimin=Pang(IP(i))-(Pang(IP(i))-Pang(IP(i-1)))*ddel/(Pdel(IP(i))-Pdel(IP(i-1)));
    phimax=Pang(IP(i))+(Pang(IP(i+1))-Pang(IP(i)))*ddel/(Pdel(IP(i+1))-Pdel(IP(i)));
    dphi=abs(phimax-phimin);
    lati=ddel;
    loni=0;
    latr=0;
    lonr=Pang(IP(i));
    theta1=heading_sphere(lati,loni,latr,lonr);
    lati=-ddel;
    theta2=heading_sphere(lati,loni,latr,lonr);
    dtheta=abs(theta2-theta1)
    dOm=dphi*pi/180*dtheta*pi/180*sin(Pang(IP(i))*pi/180);
    dOmega(i)=dOm;
    end
end

QP=2000;
QR=400;
RE=4E7./(2*pi);
V=1500;

I=find(freqs > 0.2);
i=I(1);
Fp3D_df=1E-16*betas^5*rhos^2/(2*pi*2*2*pi*freqs(i)*kdisp(i,1,9).^2);
SDF=Fp3D_df.*(2*pi*2*2*pi*freqs(i)*kdisp(i,1,9).^2)./(betas^5*rhos^2);
dA=1E10;

var_delta_P=4.*Fp3D_df/(rhos^2*betas^4)*PSdir2(i,31)*dOmega.*exp(-2*2*pi*freqs(i).*Ptau(IP)./QP);
var_delta_R=SDF./(RE.*sin(Pdel(IP)*pi/180)).*exp(-2*2*pi*freqs(i).*RE.*Pdel(IP)*pi/180./(V.*QR)).*dA;

semilogy(Pdel(IP),1E6.*sqrt(var_delta_P),'k--',Pdel(IP),1E6.*sqrt(var_delta_R),'r-','LineWidth',2)
axis([10 60 0.0001 0.1])
set(gca,'LineWidth',2,'FontSize',16,'Position',[0.1 0.165312 0.86 0.756978]);

legend('direct P phase','Rayleigh, Q=400, V=1500 m/s')
 title('Ground motion at a distance from a Fp3D*df = 4.2E4 hPa^2 m^2 source over a 100 km-side square')
  xlabel('\Delta (degrees)')
ylabel('RMSVGD (microns)')