clear
hold off
   
              %  krakcd -  normal mode code based on KRAKEN  (dev)
              %           solves eigenproblem using C routines
              %           J. Ianniello - 20 June 94



                      % INPUT PARAMETERS

 
 Nmodesmx = 44;   %  max number of modes to use and/or save
   Nmodes = Nmodesmx;       %  number of propagating modes
 
        
	fl=10;         %  lower frequency limit
        fu=10;         %  upper frequency limit
              
        df=10;           %  frequency increment
          nfreq = fix((fu - fl)/df) + 1;
        
  printkm = 1;          %  =1 if want to print wavenumbers
  plotproploss = 0;     %  = 1 if want to plot proploss
  plotlossallzs = 0;    %  = 1 if only want proploss for all receiver depths
  saveeigs = 0;         %  = 1 if want to save eigvects and vals; change name
  plotmodefuncts = 0;   %  = 1 if want to plot modefunctions
  savall = 0;           %  = 1 if want to save evals and evects for all freqs
  transfercalc=0;       %  = 1 if want to calculate transfer function

    
    if plotproploss == 1

         %  parameters for proploss

         zsm= 500;         %  source depth      
         zrm= 2500;         %  receiver depth
         r=100:25:50000;  %  ranges 

    end   %  end plotproploss set up

  
   if transfercalc==1      %  parameters for transfer function calc
         r=2000;  zsm = 31;  zrm = 38;        
   end        %  end transfercalc set up


         
           % *******   INPUT ENVIRONMENTAL MODEL    ***********
  
  
                  %  read in svp and set up environmental data

   plottitle = 'isovelocity'; % 'Gulf of Mexico';
       
   L = 500;     %  number of layers in water column;
                %  interface at L+1
                %  L should be L > 10*D/lam = 10*D*f/c
                           
 
   load svpiso.1; % svpgom1.1;  % svpgom2mod.1    % load svp and depths in water
         svp= svpiso; % svpgom1;  % svpgom2mod;      
   
              
   rhow = 1;           % water density
   rhos = 2; % 1.85;   % sediment density
   rhor = 2; %         % rock density

   csed =  2000; % 1560;  % sediment sound velocity
   crock = 2000;  %      % rock sound speed

   Dsed =  1000;  % 20;  %  depth of sediment layer
   Drock =  0;        %  depth of rock layer; = 0 if no rock layer

   attsed = .0025/(20*log10(exp(1)));  
                % compressional wave attenuation in sediment 
                % dB/m converted to nepers/m

   attrock = attsed;  % compressional wave attenuation in rocks

   sigs = .1;         %  rms roughness of surface
   sigb = .1;         %  rms roughness of bottom
       sigssq = sigs*sigs;
       sigbsq = sigb*sigb;

        rho1=rhow;
        rho2=rhos;
        rho3=rhor;
        densfact   = sqrt(2/(1+rho1/rho2));
            dfsq = densfact*densfact;
        densfactsr = sqrt(2/(1+rho2/rho3)); 
            dfsrsq = densfactsr*densfactsr;
          
            % put data in depth and c(z) arrays

              ndep=length(svp);
                             
                z=zeros(ndep,1);
                cw=zeros(ndep,1);
              z(:)=svp(:,1);    % z(1)=0; z(ndep)=D
              cw(:)=svp(:,2);    % svp in water
              D=z(ndep);      % D is depth of water column
                       
             dd = [z' D+1 D+Dsed]/D;  % augment depths with sed depth  
                 ldd=length(dd);
                 co=max(cw);

            ctot=[cw' csed csed csed]';  % add sed speed just below boundary
             
            ddtot=[dd(1:length(cw)) 1.0001 dd(length(cw)+1:ldd)];
                                      %  add a depth just below boundary

                    %  add rock layer
    if Drock > 0

      ddtot=[D*ddtot D+Dsed+.001 D+Dsed+1 D+Dsed+Drock]/D; % add rock depths
             ldd = length(dd);    
      ctot=[ctot' crock crock crock]'; % add rock speed below boundary

    end


            %   plot non-interpolated profile 

               plot(ctot,-ddtot*D);
                 title(plottitle);
                 xlabel('Sound Velocity (m/sec)');
                 ylabel('depth (m)');

               disp('pausing for SVP');
                 pause


          %  interpolate c

               d=0:1/L:ddtot(length(ddtot));  % normalized depth at int values                    
               c=ctot/co;                % normalized svp
               cint=zeros(1:length(d));
               cint(1)=c(1);
               cint(length(d))=c(length(c));

         for j=2:length(d)-1

  	   ind=find(d(j)>ddtot);
           indlow=ind(length(ind));
           dlow=ddtot(indlow);
           clow=c(indlow);
           cint(j)=clow+(d(j)-dlow)*(c(indlow+1)-clow)/(ddtot(indlow+1)-dlow);

        end
                   %  plot interpolated profile

             plot(c,-ddtot*D,cint,-d*D);
                title(plottitle)
                xlabel('Sound Velocity (m/sec)');
                ylabel('depth (m)');
 
             disp('pausing for profile and interpolated profile')
                  pause

            N=length(d)-1;    %  total number of layers
            h=1/L;
                   if Drock > 0
                       Lr = L + Drock*L/D;  %  location of rock layer
                   else
                       Lr = N;
                   end

            c=cint';
            d=d';
     
 
        if savall==1
           normphisv = zeros(N+1,Nmodesmx*nfreq);
           kmsv = zeros(Nmodesmx,nfreq);
        end


            %  *******  START FREQUENCY LOOP ************


    csq=c.*c;
    kcoef=(2*pi*D*h/co)^2;   % coefficient for k
      
   disp('  ')
   disp('  *** modesblpl running ***')
   disp('  ')
   disp('       run run run')


fc=1;
for f=fl:df:fu
                   %  f         %  print f
  ksq=kcoef*f*f;   

  kmin=2*pi*f/(co*c(N+1));  %  make sure no modes reach rigid bottom
  
 
        %   ************   SET UP A MATRIX    ******************

 
          a=-2*ones(N,1)+ksq*ones(N,1)./csq(2:N+1);
          lowdiag = ones(N,1);
          updiag = lowdiag;
       
      %  modify for interface and generalized eigenvalue problem
   
  a(L)= (-1-rho1/rho2+.5*ksq*(1/csq(L+1)+rho1/(rho2*csq(L+2))))*dfsq;  
 
             updiag(L+1)  = rho1*densfact/rho2;   % changes for sed layer
             lowdiag(L+1) = densfact;
             updiag(L)    = densfact;
             lowdiag(L)   = densfact;
     
          if Drock > 0
  a(Lr)=(-1-rho2/rho3+.5*ksq*(1/csq(Lr+1)+rho2/(rho3*csq(Lr+2))))*dfsrsq;

             updiag(Lr+1)  = rho2*densfactsr/rho3; % changes for rock layer
             lowdiag(Lr+1) = densfactsr;
             updiag(Lr)    = densfactsr;
             lowdiag(Lr)   = densfactsr;
          end

        lowdiag(N) = 2;           %  rigid botom bc
       
               %   find eigenvalues and eigenvectors

      diagsq = [updiag.*lowdiag];
      evals = trival(0,max(a)+2,a,diagsq,Nmodesmx);

      evects = trivec(evals,a,diagsq);   % evects of quasi-symm matrix
     
       
              %  convert to evects of asymm matrix

      ddiag=1;
      for ll = L+1:N
        ddiag = ddiag*sqrt(lowdiag(ll)/updiag(ll));
        evects(ll,:) = evects(ll,:)*ddiag;
      end
 
      evects(L,:)=evects(L,:)*densfact;   %  make generalized eigenvector
     
      if Drock > 0
         evects(Lr,:)=evects(Lr,:)*densfactsr; 
      end

               %   order evects and evals

           
    munorm=evals/(h*h*D*D);
    km=sqrt(munorm);

     if Nmodesmx > length(find(km>=kmin));
         Nmodes = length(find(km>=kmin));
     else 
         Nmodes = Nmodesmx;
     end

   
             % ******  PLOT NUMERICAL SOLUTIONS  *********

  for i=1:Nmodes
      depth(:,i)=d;  
  end

  ep=[zeros(1,Nmodes); evects(:,1:Nmodes)];  % put zeros in top row of evects

               
            
              % *********   PRINT HOR WAVENUMBERS  ***********

      if printkm ==1

            disp('horizontal wavenumber')
                       km'
                    
           disp('number of propagating modes')
                    num_modes=Nmodes
       end


            %  ********  NORMALIZE EIGENVECTORS ******
          
      
        watersum=sum(evects(1:L-1,:).^2);
        watersum=watersum+.5*evects(L,:).^2;
        watersum=D*watersum/L;
   
          sedsum=sum(evects(L+1:Lr-1,:).^2);
          sedsum=sedsum+.5*(evects(L,:).^2 + evects(Lr,:).^2);
          sedsum=D*sedsum/L;

        if Drock > 0
          rocksum=sum(evects(Lr+1:N-1,:).^2);
          rocksum=rocksum+.5*(evects(Lr,:).^2 + evects(N,:).^2);
          rocksum = D*rocksum/L;
          normv=sqrt(watersum/rho1+sedsum/rho2+rocksum/rho3);
        else
          normv=sqrt(watersum/rho1+sedsum/rho2);
        end
              
       normphi=zeros(N+1,Nmodes);

       for modes=1:Nmodes
         normphi(:,modes)=ep(:,modes)/normv(modes);
       end

  if plotmodefuncts==1
        hold off

       plot(normphi,-D*depth,'-r',[0 0],[-1.2 0],'-w')
       hold 
       plot([0 0],[-120 0],'w')
       plot([-.2 .2],[-31 -31],'--')           %  source depth
       plot([-.2 .2],[-38 -38],[-.2 .2],[-44 -44],[-.2 .2],[-52 -52], ...
                   [-.2 .2],[-57 -57],[-.2 .2],[-61 -61])
        title('normalized mode functions');
        xlabel('mode value')
        ylabel('depth')
        print1=sprintf('fo = %5.0f',f);
        text(.8,.95,print1,'sc');  


          disp('pausing for normalized mode functions')
          pause

  end   %  plot modefuncts==1



             %  ********  PLOT PROPAGATION LOSS

    
             %  ********   DO MODE SUM

         i=sqrt(-1);          
         sqk=sqrt(km);

             %  attsed = .0028*f^1.6/(8.6859*pi*pi);
                         
          gamBs = attsed*(2*pi*f./(rho2*co*c(L+2)*km));
         gamB=gamBs.*sedsum./(normv.*normv); 

        if Drock > 0
          gamBr = attrock*(2*pi*f./(rho3*co*c(Lr+2)*km));
         gamB=gamB+gamBr.*rocksum./(normv.*normv); 
        end
              gs =sqrt( (2*pi*f/(co*c(1)))^2-km.*km); 
              gss= sigssq*gs./(2*rho1*km); 
            gamSS = gss.*normphi(2,:).*normphi(2,:)/h^2;               
        
              g1 =sqrt( (2*pi*f/(co*c(L+1)))^2-km.*km); 
              g2 =sqrt( km.*km-(2*pi*f/(co*c(L+2)))^2 );
              gbs= sigbsq*g1./(2*rho1*rho2*km); 
              gbs = gbs.*(rho2*g1.^2+rho1*g2.^2);
            gamBS = gbs.*normphi(L+1,:).*normphi(L+1,:); 


    if plotproploss == 1 | transfercalc ==1 | plotlossallzs==1
 
            zs=fix(zsm*L/D +.5) + 1;
            zr=fix(zrm*L/D +.5) + 1;

      con=i*4*pi*exp(-i*pi/4)/sqrt(8*pi);
   
      p=zeros(1,length(r));     
      temp=normphi(zs,:)./sqk(1:Nmodes);
      expon=exp((i*km'-gamB')*r);
   
      p = (temp.*normphi(zr,:))*expon(1:Nmodes,:);
     
      if transfercalc == 1
          H(fc)=p;
          fc=fc+1;
      end

 if plotproploss==1
      clg
      hold off

       %  nummd=10;     %  to plot only a few modes
       % p = (temp(1:nummd).*normphi(zr,1:nummd))*expon(1:nummd,:);

         p=con*p./sqrt(r);
         axis([0 57000 -114 -60]);
        plot(r,10*log10(p.*conj(p)),r,-10*log10(r)-50)
          title(plottitle)
          xlabel('range (m)')
          ylabel('propagation loss (db)')
 
          print1=sprintf('fo = %5.0f',f);
          text(.8,.85,print1,'sc');    
          print2=sprintf('SD = %5.0f',zsm);
          text(.8,.8,print2,'sc');
          print3=sprintf('RD = %5.0f',zrm);
          text(.8,.75,print3,'sc');

disp('pausing for proploss')
  pause

  if plotlossallzs==1
   hold off
   clg
   disp('  ')
   disp('calculating prop loss for all r and z')
 
      p=zeros(L,length(r));
   
      for ll=1:L
        normphi(ll,:)=normphi(ll,:).*temp;
      end
      p=normphi(1:L,:)*expon(1:Nmodes,:);
  
      p=con*p;
      for l=1:L
        p(l,:)=p(l,:)./sqrt(r);
      end
   
      p=10*log10(p(2:L,:).*conj(p(2:L,:)));
         mxp=max(max(p));
         p=p-mxp;
         
         %  v=[-5 -10 -15 -20 -25 -30 -35 -40 -45 -50];
         %  v=[-10 -20 -30 -40 -50 -60]
            v = 10;
        contour(p,v,r,[-D:D/L:-2*D/L])
          

        title(plottitle)
        xlabel('range (m)')
        ylabel('depth (m)')

        print1=sprintf('fo = %5.0f',f);
        text(.8,.95,print1,'sc');  
        print2=sprintf('SD = %5.0f',zsm);
        text(.8,.9,print2,'sc');

  end   % plotlossallzs
 end    % plot proploss   

 end   % end if all...
 
       if savall == 1
         normphisv(:,Nmodesmx*(fc-1)+1:Nmodesmx*(fc-1)+Nmodes)=normphi;
         kmsv(1:Nmodes,fc)=km;
         fc = fc + 1;
       end
           
 end  % f loop

                   %  save eigenvalues and vectors
   if savall == 1
     normphisvsv = normphisv(1:L,:);
     save gomsvp2mod normphisvsv kmsv 
   end

if saveeigs==1
  save eigvblpl2 normphi;
  %   put in save of km
end
                  %  plot transfer function

if transfercalc==1
   freq=fl:df:fu;
   plot(freq,abs(H)/max(abs(H)));
   title('Transfer Function')
   xlabel('frequency (Hz)')
   ylabel('Magnitude of Transfer Function')
end
