 
        % bart   -  loads replicas and does bartlett matched field
        % 	processing using kraken transfer functions  
                
        %       J. Ianniello   -    12 Sept


clear
clg
hold off
 
  
       			% INPUT PARAMETERS

  
 	     zr = [38 41 52 57 61];   %  known line depths cs8a
           % zr = [121 145 165 180 195]*.305;  % cs5a
           % zr = [122 145 168 182 196]*.305;  % cs3a
           % zr = [112.5 135 156 171 183]*.305; % cs6a
             
        nlines = length(zr);
             linest =1;
             linestp = nlines;

	   
	fl=100;      % lower frequency
	fu=600;      % upper frequency
	df=1        % frequency increment
			          
        fpoints = fix((fu-fl)/df+1);   % # of frequency points

    
              		% SET UP SEARCH REGION


	 
	D = 111;                     
	dpoints=400;        % # of depth points
           L = dpoints;
        rtrue = 21500;
        dtrue = 37;
	rmin=15000;            % minimum search range
	rmax=30000;            % maximum search range
	dr=25;                  % range search increment
        r = rmin:dr:rmax;
	      rpoints = length(r)        

                 % input model info
 
 	                             
            i = sqrt(-1);  
            zri = fix(zr*L/D +.5) + 1;
               
                  

	              
       			 %   READ REPLICAS


     realdat = 1     %  note  *** change numtmp for data or sim

            S = zeros(fpoints,nlines);
      if realdat == 1      %  read real data

            R = zeros(3125,5);
                 % load Ravgn8b.mat
                    % load Rsv5a_1_8sec.mat;  % Ravgn8b 
                    %  load Rsv3a_8cl.mat
                    % load Rsv3a_10.mat
                    % load Rsv6a_1cl.mat
                      load Rsv8b_2cl.mat
            plottitle = 'CS8b - Bartlett Processing: df = 1 Hz';

          for line =1:nlines
           R(1:2000,line) = Rss(1:2000,line)-(5-line+.5);
          end

         plot(R)
         disp('pause to plot data');
          pause

                   St = fft(R);
                  % St = ifft(R);


            for fc = 1:fpoints   
                   strt = round(fl+df*(fc-1)-df/2+.001);
                   stp  = round(fl+df*(fc-1)+df/2-.001);
               if df == 1     
                 S(fc,1:nlines) = St(strt,1:nlines);
               else
    		 S(fc,1:nlines) = sum(St(strt:stp,1:nlines));
               end
            end
          else  % do simulation
       
           % load TF_D108_df1.mat
            % load TF_D113_df1.mat
            load TF_D115_R2093_df1.mat
           S = TF;
           S = conj(S);     %  must do for simulated data only
        plottitle = 'Simulated data - Bartlett Processing: df = 1 Hz';
     end            %  end realdat
 



                                 %   LOAD STORED MODEL DATA 

          disp('loading model')

            nmodes = 25;      
            nmodusemx = 25;   %  max no of modes to use

          % load svpgom1_D108_df1.mat  % D=108; L=300; df=1
          % load svpgom2m_D113_df1.mat   % D = 113 L = 300 df = 1
          % load svpgom2m_D115_df1_20md.mat  % D = 115
          % load svpgomcs3_D115_df1.mat  % D = 115
          % load svpgomcs6_D122_df1.mat
          % load svpgomcs6_D125_df1.mat
          %  load svpgomcs6_D121_moreatt_df1.mat
             load svpgomcs5_D111_df1.mat

                 size(normphisv)

                 % strtg = 1;
                   % for gct = strtg:nlines
                      % filetoread = sprintf('load Ravg8b_%1.0f',gct);
                      % eval(filetoread)
                      % R(:,gct-strtg+1)=Ravg(1:3125);
                   % end


             % COMPUTE BART AT EACH R AND D FREQUENCY BY FREQUENCY

  	%  BART = zeros(dpoints-1,rpoints);
                          

         %	 Find transfer functions for all line depths
         %       and ranges at this frequencies

         Lstp = 4;
         L3 = L/Lstp;

             T =      zeros(L3-1,rpoints);
             BART =   zeros(L3-1,rpoints);
             numtmp = zeros(L3-1,rpoints);
             dentmp = zeros(L3-1,rpoints);

          % TFF=zeros(fpoints,nlines);
          
  for fc = 1:fpoints                                      	    
     fc  
    
  	kmt = kmsv(1:nmodes,fc);
         nmoduse = length(find(kmt>0));   %  number of non-zero modes
              if nmoduse > nmodusemx
                  nmoduse = nmodusemx;
              end
        km = kmt(1:nmoduse);
        gamB = gamBsv(1:nmoduse,fc);
         sqk = sqrt(km(1:nmoduse));
         normphi = zeros(L3,nmoduse);

normphi(1:L3,:)=normphisv(2:Lstp:L-1,nmodes*(fc-1)+1:nmodes*(fc-1)+nmoduse);
        % normphi(L3-1,1:nmoduse) = normphi(L-2,1:nmoduse);

        temp2 = zeros(L3-1,nmoduse);                         
	expon = exp((i*km-gamB)*r);
           
   for lines = linest:linestp; %1:nlines
       
      temp = normphi(round(zri(lines)/Lstp),:)./sqk';
             
      for ll=1:L3-1
        temp2(ll,:)=normphi(ll,:).*temp;
      end
      
      T=temp2*expon;
                     
       numtmp = numtmp + S(fc,lines)*T;
       dentmp = dentmp + conj(T).*T;
                                             % TFF(fc,lines) = T(112);
   end  %   end nlines
 

     BART = BART + real(conj(numtmp).*numtmp./dentmp); % incoherent
             numtmp = zeros(L3-1,rpoints);
             dentmp = zeros(L3-1,rpoints);
 
  end     %  end f loop

    %   BART =  real(conj(numtmp).*numtmp./dentmp);  % coherent
             

		%  PLOT RANGE/DEPTH

 clg
 hold off

  mesh(BART)
   
    title(plottitle)
    text(.8,.2,'range (km)','sc');
    text(.05,.3,'depth (m)','sc');
         labrmn = sprintf('%6.0f',rmin);
         labrmx = sprintf('%6.0f',rmax);
         labd   = sprintf('%5.0f',D);
    text(.5,.15,labrmn,'sc')
    text(.9,.31,labrmx,'sc')
    text(.4,.16,labd,'sc')  
    text(.06,.4,'0','sc')

    lab1 = sprintf('true range = %4.1f',rtrue);
    text(.1,.9,lab1,'sc')
    lab2 = sprintf('true depth = %4.1f',dtrue);
    text(.1,.85,lab2,'sc')
     
     [bigvb,indvb] = max(BART);  % max in columns; row indices of max`s
     [bignb,indcb] = max(bigvb);    % biggest #, col index of max 
 		indrb=indvb(indcb);
 
      	depestb = D*indrb/L3
	ranestb=dr*(indcb-1)+rmin

     printr = sprintf('range peak = %4.1f',ranestb);
     text(.65,.9,printr,'sc')
     printd = sprintf('depth peak = %4.1f',depestb);
     text(.65,.85,printd,'sc')

     % save BARTcs8b2_cl_111.mat BART rmin rmax D rtrue dtrue L3 dr plottitle
       
  
