clear
hold off
 
  
              %  krakc  -  normal mode code based on KRAKEN   
              %             solves eigenproblem using C routines
               
              %           J. Ianniello - 12 Sept 94       
                     
              %            C routines trival and trivec by W.R. Bernecky
               


                      % INPUT PARAMETERS

 
 Nmodesmx = 20;   %  max number of modes to use and/or save
   Nmodes = Nmodesmx;       %  number of propagating modes
 
        
	fl=100;         %  lower frequency limit
        fu=600;         %  upper frequency limit
              
        df=3;           %  frequency increment
          nfreq = fix((fu - fl)/df) + 1;
        
  printkm = 0;          %  =1 if want to print wavenumbers
  plotmodefuncts = 0;   %  = 1 if want to plot modefunctions

  plotproploss = 0;     %  = 1 if want to plot proploss at 1 rec depth
  plotlossallzs = 0;    %  = 1 if want proploss for all receiver depths
  transfercalc= 0;       %  = 1 if want to calculate transfer function
    plotimp = 0;

  savall = 1;           %  = 1 if want to save evals and evects for all freqs
  
    
    if plotproploss == 1  |  plotlossallzs == 1

         %  parameters for proploss

         zsm=  40;         %  source depth      
         zrm=  36;         %  receiver depth; put in proper range for allzs
         r = 2500:10:20000;  %  ranges 

    end   %  end plotproploss set up

  
   if transfercalc==1      %  parameters for transfer function calc
          
              r = 1852*5.6;  zsm = 80*.305; zrm = [140 143 146 149 154]*.305;
               r=9000;
           datatitle = 'V3S1';  

          nlines = length(zrm);
          H = zeros(nlines,nfreq);
   end        %  end transfercalc set up


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

   plottitle = 'Propagation Loss for GOM';
       
   L = 400;     %  number of layers in water column;
                %  interface at L+1
                %  L should be L > 10*D/lam = 10*D*f/c
                           
 
   load  gom2b94.dat;   %svpgom2_111.1;  %  load svp and depths in water
         svp = gom2b94;  %svpgom2_111;  %                       
   rhow =  1;           % water density
   rhos =  1.9;         % sediment density
   rhor =  2.2;         % rock density

   csed =  1505; %1560;  % sediment sound velocity
   sedratio = 1.1 ;  %  ratio of csed(0) to cw(D); = 0 if dont want to use
   sedprof = 1;          %  = 1 for constant or linear; =2 for Hamilton
      sgd = 0;          % gradient of sediment sound speed (m/sec)/m
   crock = 3000  % rock sound speed

   Dsed = 2.5;  %3.5;  %  depth of sediment layer
   Drock = 10;        %  depth of rock layer; = 0 if no rock layer

   attsednpl =.56/(20*log10(exp(1)));   % compressional wave attenuation 
                          % sediment  dB/lam converted to nepers/lam

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

        rho1=rhow;
        rho2=rhos;
        rho3=rhor;

        dfsq   =  2/(1+rho1/rho2);          
        dfsrsq =  2/(1+rho2/rho3); 
                       
            % 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
              
          if sedprof == 1         
             dd = [z' D+1 D+Dsed]/D;  % augment depths with sed depth  
                 ldd=length(dd);
                 co=max(cw);

              if sedratio ~= 0
                 csed = sedratio*cw(ndep);  % if want sed speed from ratio
              end
 
            ddtot=[dd(1:length(cw)) 1.0001 dd(length(cw)+1:ldd)];
                                      %  add a depth just below boundary
            ctot=[cw' csed csed+sgd csed+sgd*Dsed]';  
                           % add sed speed just below boundary

          else
            Dsed = Dsed + .05;
            ddtot = [z' D+.05:1:D+Dsed]/D;  % augment depths with sed depth   
            co=max(cw);
            K = sedratio*cw(ndep)/(.05^0.015);
            ctot=[cw' K*[.05:1:Dsed].^.015]';  
         end              
             
            
                    %  add rock layer
    if Drock > 0

      ddtot=[D*ddtot D+Dsed+.001 D+Dsed+1 D+Dsed+Drock]/D; % add rock depths   
      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 = fix(L+Dsed*L/D);  %  location of rock layer
                   else
                       Lr = N;
                   end

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


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


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

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

  kmin = 2*pi*f/(co*c(Lr+2));  %  make sure no modes reach rigid bottom

  attsed = attsednpl*f/csed;   % compressional wave attenuation in sediment 
                               % nep/lam converted to nepers/m
   attrock = .5*attsed;  % compressional wave attenuation in rocks

  
 
        %   ************   SET UP DIAGONALS OF A   ******************

 
          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*dfsq/rho2;   % changes for sed layer
             lowdiag(L)   = dfsq;
      
          if Drock > 0
  a(Lr)=(-1-rho2/rho3+.5*ksq*(1/csq(Lr+1)+rho2/(rho3*csq(Lr+2))))*dfsrsq;

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

        lowdiag(N) = 2;           %  rigid botom bc
       

               %   find eigenvalues and eigenvectors

      diagsq = [updiag.*lowdiag];
      evals = trival(0,max(a)+2,a,diagsq,Nmodesmx);
            
  
               %   scale evals
           
    km = sqrt(evals/(h*h*D*D));
    
     if Nmodesmx > length(find(km>=kmin));
         Nmodes = length(find(km>=kmin));
     else 
         Nmodes = Nmodesmx;
     end

     km = km(1:Nmodes);
     evals = evals(1:Nmodes);

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

      ddiag=1;
      for ll = L:N
        ddiag = ddiag*sqrt(lowdiag(ll)/updiag(ll));
        evects(ll,:) = evects(ll,:)*ddiag;
      end

            
   % *********   PRINT HORIZONTAL WAVENUMBERS AND PHASE SPEED ***********

      if printkm ==1

            disp('    horizontal wavenumber')
              disp(' ')
            disp('       Kr                 phase speed')
                      [ km'  2*pi*f./km']
             
           disp('number of propagating modes')
                    num_modes=Nmodes
       end


            %  ********  NORMALIZE EIGENVECTORS ******
          
                   %   get integrals
      
        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)= [0; evects(:,modes)]/normv(modes);
       end



                 %   ***  PLOT MODEFUNCTIONS  *****


  if plotmodefuncts==1
        hold off
        clg
                     
        for i=1:Nmodes
          depth(:,i)=d;  
         end

       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

          lossum = zeros(1,Nmodes);     
         for m=1:Nmodes
           lossum(m)=sum((evects(L+1:Lr-1,m).^2)./c(L+2:Lr));
           lossum(m)=lossum(m)+.5*( (evects(L,m).^2)/c(L+2) ...
                         +(evects(Lr,m).^2)/c(Lr+1) );
         end    
          lossum=D*lossum/L;
                         
           gamBs = attsed*(2*pi*f./(rho2*co*km));
         gamB  = gamBs.*lossum./(normv.*normv); 
 
        if Drock > 0
           gamBr = attrock*(2*pi*f./(rho3*co*c(Lr+2)*km));
         gamB=gamB+gamBr.*rocksum./(normv.*normv); 
        end
              gs =real(sqrt( (2*pi*f/(co*c(1)))^2-km.*km)); 
              gss= sigssq*gs./(2*rho1*km); 
            gamSS = gss.*normphi(2,:).*normphi(2,:)*L^2/D^2;               
        
              g1 =real(sqrt( (2*pi*f/(co*c(L+1)))^2-km.*km)); 
              g2 =real(sqrt( km.*km-(2*pi*f/(co*c(L+2)))^2));
              gbs= sigbsq*g1./(2*rho1*rho2^2*km); 
              gbs = gbs.*(rho2^2*g1.^2+rho1^2*g2.^2);
            gamBS = gbs.*normphi(L+1,:).*normphi(L+1,:); 

            gamB = gamB + gamBS + gamSS;

     i=sqrt(-1); 

    if plotproploss == 1 | transfercalc ==1 | plotlossallzs==1

            sqk=sqrt(km);
            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);
         
      temp=normphi(zs,:)./sqk(1:Nmodes);
      expon=exp((i*km'-gamB')*r);
       
      if transfercalc == 1
         for line = 1:nlines
          H(line,fc)=con*(temp.*normphi(zr(line),:))*expon(1:Nmodes,:)/sqrt(r);
         end
          fc=fc+1;
      end

 if plotproploss==1
      clg
      hold off
           disp('  ')
           disp('calculating prop loss vs r')

          p=zeros(1,length(r)); 
          p = (temp.*normphi(zr,:))*expon(1:Nmodes,:);

       %  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 70000 -140 -50]);
           axis([0 20000 -90 -60]);
        plot(r,10*log10(p.*conj(p)));  % ,r,-20*log10(r))
          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

  end    % plot proploss  

  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,:)));
                  
          % v=[-5 -10 -15 -20 -25 -30 -35 -40 -45 -50];
          % v=[-10 -20 -30 -40 -50 -60]
          v = [-45 -50 -55 -60 -65 -70 ];
             
        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   % end if all...
 

         %   save modefunctions and hor wavenumbers

       if savall == 1
         normphisv(:,Nmodesmx*(fc-1)+1:Nmodesmx*(fc-1)+Nmodes)=normphi;
         kmsv(1:Nmodes,fc)=km';
         gamBsv(1:Nmodes,fc)=gamB';
         fc = fc + 1;
       end
           
 end  % *******  END f LOOP  ************

 

                   %  store modefunctions and hor wavenumbers
     if savall == 1
         disp('double pause at savall')
         pause; pause;  disp('saving model info');
       % normphisvsv = normphisv(1:L,:);
        save svp_D94_df3 normphisv kmsv gamBsv
       disp('saved normphi')
     end

     
                  %  plot transfer function

  if transfercalc==1
     lll=2;
     freq=fl:df:fu;
     plot(freq,abs(H(lll,:)));
     title('Magnitude of Transfer Function')
     xlabel('frequency (Hz)')
     ylabel('Magnitude')
      print1=sprintf('range = %5.0f',r);
      text(.73,.89,print1,'sc');  
      print1=sprintf('source depth = %5.0f',zsm);
      text(.73,.84,print1,'sc'); 
      print1=sprintf('receiver depth = %5.0f',zrm(lll));
      text(.73,.79,print1,'sc'); 

     disp('pausing for transfer function')
     pause


   if plotimp == 1

          
       nplt=1500; %1400; %512;
     rcop = zeros(nplt,nlines);
         
           load Rsv_v3s1_1.mat

       rsv = zeros(length(Rss),nlines);
       for line = 1:nlines
         rsv(:,line) = Rss(:,line) - (5-line+.5);
       end

              lll=1;
       fs = 3125;         %  assumes df = 1 Hz
       Nf = 3125;  freq=fl:df:fu;
    mrsv=max(abs(rsv));
    rcop = zeros(nplt,1);
       nst =75 ;   %250; %160; % 334;
      rcop(nst:nplt)=rsv(1:nplt-nst+1,lll)/mrsv(lll);
      
       % rcop(1:nplt)=rsv(nst:nplt+nst-1,lll)/mrsv(lll);

       HH = zeros(3125,1);
        
        HH(fl+1:fu+1) = H(lll,1:nfreq).*exp(-i*2*pi*[fl:fu]*(r/co-.05));
        HH(Nf:-1:(Nf+1)/2+1)=conj(HH(2:(Nf+1)/2));   
   
       impres = real(fft(HH));
       max(abs(impres))
       impres=impres/max(abs(impres));
      % plot( [1:nplt]/fs,impres(1:nplt),[1:nplt]/fs,rcop(1:nplt));
        % plot( [1:nplt]/fs,impres(1:nplt));

    clg
    hold off
    subplot(211)
         plot( [1:nplt]/fs,impres(1:nplt));
          title('simulation')
          print1=sprintf('range = %5.0f',r);
          text(.7,.91,print1,'sc');  
          print1=sprintf('source depth =%5.0f',zsm);
          text(.7,.87,print1,'sc'); 
          print1=sprintf('receiver depth =%5.0f',zrm(lll));
          text(.7,.83,print1,'sc'); 
          print1=sprintf('water depth =%5.0f',D);
          text(.7,.67,print1,'sc'); 
          print1=sprintf('sedratio =%5.2f',sedratio);
          text(.7,.63,print1,'sc');

    subplot(212)
         plot([1:nplt]/fs,rcop(1:nplt));
           title(datatitle)
           print1=sprintf('line = %5.0f',lll);
           text(.8,.41,print1,'sc');
           print1=sprintf('time = %5.0f',corrtimeh);
           text(.8,.38,print1,'sc');
           print1=sprintf('                %5.0f',corrtimem);
           text(.8,.381,print1,'sc');


    end   %  end plot imp


  end    %   end transfer calc
 
