function TLm=TLm_calc(fft_root,out_root,Tlims,dT,ph_lims,flims,ii_phavg,SL,SL_parms,iiAN_corr,AN_freqs, ...
    geom_file)
%
% Computes measured TL from FFT file and SL, measured SL, or measured AN.
%
% Call TLm=TLm_calc(fft_root,out_root,Tlims,dT,ph_lims,flims,ii_phavg,SL,SL_parms,iiAN_corr,AN_freq,geom_file);
%
% fft_root	= Input FFT file root;
% out_root	= Output file root for storing TLm structure ([]=do not save);
% Tlims		= Time interval [Tmin,Tmax] to process (in HOD) ([]=all);
% dT		= Time interval to process at one time (0 or [] for auto calc);
% ph_lims	= List of receivers to process;
% flims		= Frequency interval(s) [fmin1,fmax1; fmin2,fmax2; ...] to process ([]=all);
% ii_phavg	= 1 to average over receivers, 0 to output receiver separately;
% SL		= Source level (single number or two columns of [t_HOD,SL])
%			  Set SL=0 for processing reference phone data for source level;
%			  Set SL=-1 for simulated TL FFT file;
%			  Set SL=-999 for to compute ambient noise in dB//uPa^2/Hz instead of TL;
%             SL = structure containing SL(t) from previous run of TLm_calc on reference hydrophone;
% SL_parms  = [SL_dyn,T_filt], where
%             SL_dyn    = dB drop to indicate source off and NaN fill should be done
%             T_filt    = Period cutoff for low-pass filter of SL data;
% iiAN_corr = Flag to correct for background AN (assumes AN lies between SL lines; uses median to find AN);
% AN_freqs  = Frequency band over which to estimate AN;
% geom_file = Array geometry file to use to fill TLm.scales2 with array depths ([] ==> do not fill);
%
% Examples:
%	TLm=TLm_calc('swami_tl','swami_tl',[6,10,.5],[],[225,243],185.8,[6,2.5])
%	TLm=TLm_calc('swami_tl_ev74','swami_tl_ev74',[16.5,20,.5],[],[225,243],185.8,[6,2.5])
%
if nargin<8, SL=0; end
if nargin<9 | isempty(SL_parms), SL_parms=[9,10]; end
if nargin<10 | isempty(iiAN_corr), iiAN_corr=0; end
if nargin<11, AN_freqs=[]; end
%
SL_dyn=SL_parms(1); T_filt=SL_parms(2);
%
if isempty(out_root), out_root=fft_root; end
% Read FFT header:
FFT=fft_read1_open(fft_root);
fclose(FFT.fft_fid);
if isempty(ph_lims), iiph=[1:FFT.nrec]'; else, iiph=adi_lim_eval(ph_lims); end
% Check if geometry file valid for FFT file:
ii_geom=0;
if ~isempty(geom_file),
	[xyz_all,cavg,nph_all]=arr_geom(geom_file);
    if any(iiph>FFT.nrec),
        error(['# phones in FFT file > # phones in array geometry file: ', ...
                int2str([FFT.nrec,nph_all]),'.  Aborting ...']);
    elseif FFT.nrec<nph_all,
        disp(['NOTE: # phones in FFT file < # phones in array geometry file: ', ...
                int2str([FFT.nrec,nph_all]),'.  Assuming first FFT.nrec phones ...']);
    end
    ii_geom=1;
end
%
if isempty(flims), flims=[FFT.f1,FFT.f2]; end
fmin=min(flims(:)); fmax=max(flims(:)); 
%
%fax=(FFT.imap(:,1)-1)*FFT.df;
%jfreq=find(fax>=fmin-FFT.df/2 & fax<=fmax+FFT.df/2);
%nfreq=length(jfreq); fax_get=fax(jfreq); 
%
[nband,nc]=size(flims);
if nc~=2, disp(flims); error('User error: frequency limits must be nband by 2 array'); end
%
if isempty(Tlims), Tlims=FFT.hodax([1,end]); end
[nr,nc]=size(Tlims);
if nc==1,
	if ~all(round(Tlims)==Tlims), 
		error('Column vectors for Tlims must be integer indices');
	end
	Tmin=FFT.hodax(min(Tlims)); Tmax=FFT.hodax(max(Tlims));
	if Tmin==Tmax, Tmax=Tmin+1; end
else	
	Tmin=Tlims(1); Tmax=Tlims(2);
end
iit=find(FFT.hodax>=Tmin & FFT.hodax<=Tmax); iit=iit(:); ntime=length(iit); 
t_ax=FFT.hodax(iit); t_ax=t_ax(:);
%
T1=Tmin;
t_hour=[]; RL=[]; nq=0;
if length(iiph)==1 | ii_phavg, nrec=1; else; nrec=length(iiph); end
%
f0=mean(flims,2); BW=flims(:,2) - flims(:,1);
% If FFT file is simulated TL (SL=-1), correct for df factor in calcs:
iiavg=2; % Default is to integrate over frequency for TL and SL calcs;
dfnorm=1;
suf='_TLm'; nfilt=1;
if isstruct(SL),
    if SL.dims(2)~=1, disp('WARNING: Unexpected number of receivers in SL structure: '); disp(SL); end
    iiSL=3; 
    Jday_SL=str2num(SL.labels{3}(end-2:end));
    t_ax_SL=(Jday_SL - FFT.jday0)*24 + SL.scales{3}(:); 
    SL_ax=zeros(SL.dims(3),nband);
    for jband=1:nband,
        [dummy,kband]=min(SL.scales{1}-f0(jband));
        SL_ax(:,jband)=reshape(SL.data(kband,1,:,end),[SL.dims(3),1]);
    end
    if t_ax_SL(1)>t_ax(end) | t_ax_SL(end)<t_ax(1),
        disp(t_ax([1,end])'); disp(t_ax_SL([1,end])');
        error('No overlap between FFT times and source level SL times.');
    end
elseif length(SL)==1,
    if SL==-1,
		disp(['For SL=-1, TLm_calc averages over modeled spectra rather than ']);
		disp(['integrates spectral levels using df'])
		iiavg=1; dfnorm=0;   % Average over frequency since modeled spectra do not depend on df
        iiSL=-1;
    elseif SL==-999,
		disp(['For SL=-999, TLm_calc computes AN in dB//uPa^2/Hz.']);
        iiavg=1; suf='_AN';   % Compute ambient noise by averaging dB//uPa^2/Hz over frequency;
        iiSL=-999;
    elseif SL==0,
        suf='_SL'; iiSL=0; nfilt=2;
    else
        iiSL=1;
    end
else
	[nr,nc]=size(SL);
	if nc~=2, 
		disp(['Illegal SL(t) specification: ', ...
		num2str(SL(:)'),' Must have two columns: t_hr,SL.']); 
		return 
	end
	if SL(1,1)>Tmin | SL(nr,1)<Tmax,
		disp('First, last times in [t_HOD,SL] array do not cover time interval requested');
		disp(['SL time interval: ',sprintf('%2.2f,2.2f',SL([1,nr],1)), ...
                '; Tlims = ',sprintf('%2.2f,2.2f',Tlims)]);
	end
    iiSL=2; t_ax_SL=SL(:,1); SL_ax=SL(:,2);
end		
% Do not integrate over f for modeled spectra;
% Average over receivers, integrate (dfnorm=1) or average (dfnorm=0) over f to get uPa^2:
kdim=1;
if ii_phavg, kdim=[kdim 2]; iiavg=[iiavg 1]; end        % Average over phone also;
if iiAN_corr==1,
    if ii_phavg, iiavg2=[4 2]; else, iiavg2=[4]; end
    BW_fac=10*log10(BW);
elseif iiAN_corr==2,
    if size(AN_freqs,1)~=nband, 
        error('USER ERROR: AN_freqs must be same size as TL_freqs for AN correction.'); 
    end
    iiavg2=1; iiavg2=4,
    if ii_phavg, iiavg2=[iiavg2 1]; end
    BW_fac=10*log10(BW);
%    AN_freqs=flims + repmat(1.5*BW,[1,2]);
%    iibad=find(AN_freqs(:,2)>FFT.f2); nbad=length(iibad);
%    if ~isempty(iibad),
%        AN_freqs(iibad,:)=flims(iibad,:) - repmat(1.5*BW,[1,2]);
%        iibad=find(AN_freqs(:,1)<FFT.f1); nbad=length(iibad);
%        if ~isempty(iibad),
%            AN_freqs(iibad,1)=FFT.f1;
%        end
%    end
%    BW_fac=10*log10(AN_freqs(:,2)-AN_freqs(:,1));
    disp('AN Frequency limits: '); disp(AN_freqs);
end
%
iihr=1; iinorm=0; parms=[dfnorm, iihr, iinorm, 1, 0];
%
if isempty(out_root), out_root=fft_root; end
out_file=[out_root,suf];
TLm.title=out_root;
TLm.fft_root=fft_root;
TLm.ndim=4; TLm.dims=[nband nrec ntime nfilt];
TLm.labels={'Center Frequency - Hz', 'Receiver Number', ['Time - Hour of Day ',int2str(FFT.jday0)], ...
        'Raw/Filtered'};
TLm.scales={f0(:)', iiph(:)', t_ax, [0:nfilt-1]};
TLm.data=NaN*ones(nband,nrec,ntime,nfilt);
%
if isempty(dT) | dT==0,
    word_max=5000000;
	word_seq=length(iiph)*FFT.nwords; dT_seq=FFT.T_step/3600;
	n_seq=max(2,word_max/word_seq); dT=n_seq*dT_seq;
end
if FFT.nseq>1, delT=diff(FFT.hodax([1,2])); else, delT=100; end
iis=max(1,round(dT/delT));
iidiag=0;
for jtime1=1:iis:ntime,
    ii=[jtime1:min(ntime,jtime1+iis-1)];
    iiseq=iit(ii);
% Read FFT portion:
    for jband=1:nband,
    	[LG,FFT]=lgram_3d(fft_root,iiseq,iiph,flims(jband,:),parms,0);
    	if isempty(LG.scales{3}), break; end
    	SP=data_avg(LG,kdim,[],iiavg,0,0);
        if iiAN_corr==1,
            MED=data_avg(LG,kdim,[],iiavg2,0,0);         % Compute Median level
            SP.data=10*log10(10.^(SP.data/10) - 10.^((MED.data + BW_fac(jband))/10));
        elseif iiAN_corr==2,
        	[LG2,FFT]=lgram_3d(fft_root,iiseq,iiph,AN_freqs(jband,:),parms,0);
            AVG=data_avg(LG2,kdim,[],iiavg2,0,0);         % Compute band-integrated noise level
            data_corr=max(10.^(SP.data/10) - 10.^((AVG.data + BW_fac(jband))/10),1e-30);
            if iidiag, SP.data(1:10) - 10*log10(data_corr(1:10)), end
            SP.data=10*log10(data_corr);
        end
        TLm.data(jband,1:nrec,ii,1)=SP.data;
    end
	disp(['Finished time ',sprintf('%2.2f',LG.scales{3}(end)),' of ', ...
            sprintf('[%2.2f,%2.2f]',[Tmin,Tmax])]);
    pause(.001);
end
%
if iiSL==2,
	SL_t=interp1(t_ax_SL,SL_ax,t_ax(:));
	TLm.data(:,:,:,1)=TLm.data(:,:,:,1) - repmat(SL_t,[nband,nrec,1]); TLm.dlab='TL - dB';
elseif iiSL==3,
    for jband=1:nband,
		SL_t=interp1(t_ax_SL,SL_ax(:,jband),t_ax);
        for jrec=1:nrec,
    	    TLm.data(jband,jrec,:,1)=TLm.data(jband,jrec,:,1) - reshape(SL_t(:,jband),[1,1,ntime]); 
        end
    end
    TLm.dlab='TL - dB';
elseif iiSL==-999,
	SL_t=SL; TLm.dlab='AN - dB//\muPa^2/Hz';
elseif iiSL==-1,
    SL_t=0; TLm.dlab='Simulated TL - dB';
elseif iiSL==0,
	SL_t=SL; TLm.dlab='SL - dB//1uPa@1m';
    dT=diff(TLm.scales{3}(1:2))*60;     % Sample period in minutes;
    f0_filt=1/T_filt; fs2_filt=1/dT/2; Wn=min(1,f0_filt/fs2_filt); 
    for jband=1:nband,
        for jrec=1:nrec,
            SL_max=max(TLm.data(jband,jrec,:,1));
            ii_NaN=find(TLm.data(jband,jrec,:,1)<SL_max - SL_dyn);
            TLm.data(jband,jrec,ii_NaN,1)=NaN;
            ii_diff=find(diff(ii_NaN)>1); ndiff=length(ii_diff);
            if isempty(ii_NaN),
                n_int=1; ii_int=[1,ntime];
            else
                n_int=0; ii_int=[];
                if ii_NaN(1)>1, n_int=n_int + 1; ii_int=[ii_int; 1,ii_NaN(1)-1]; end
                for j_int=1:ndiff,
                    jdiff=ii_diff(j_int);
                    n_int=n_int + 1; ii_int=[ii_int; ii_NaN(jdiff)+1, ii_NaN(jdiff+1)-1];
                end
                if ii_NaN(end)<ntime, n_int=n_int + 1; ii_int=[ii_int; ii_NaN(end)+1,ntime]; end
            end
%%            disp(ii_int),
            for j_int=1:n_int,
                ii=[ii_int(j_int,1):ii_int(j_int,2)]'; nii=length(ii);
                N=floor(min((nii-1)/3,10)); 
                if N==0,
                    TLm.data(jband,jrec,ii,2)=TLm.data(jband,jrec,ii,1);
                else
                    [b,a]=butter(N,Wn);
                    data=squeeze(TLm.data(jband,jrec,ii,1));
                    TLm.data(jband,jrec,ii,2)=filtfilt(b,a,data);
                end
%				delf_start=mean(delf(1:4,:)); delf_end=mean(delf(end-3:end,:));
%				delf2=[delf_start; delf; delf_end]; delf3=filtfilt(b,a,delf2);
%				delf=delf3(2:end-1,:);
            end
        end
    end
elseif iiSL==1,
	SL_t=SL; TLm.data(:,:,:,1)=TLm.data(:,:,:,1) - SL_t; TLm.dlab='TL - dB';
end
%
cmax=ceil(max(TLm.data(:))/5)*5;
cmin=max(cmax-60,floor(min(TLm.data(:))/5)*5);
TLm.clim=[cmin,cmax];
TLm.iidB=1;
%
%TLm.Tlims=[Tmin,Tmax]; 
%TLm.dT=dT;
TLm.flims=flims;
TLm.ii_phavg=ii_phavg;
%TLm.jrec=iiph;
%TLm.nrec=size(RL,1);
TLm.SL=SL_t;
if iiSL==3, TLm.SL_root=SL.title; end
TLm.jday=FFT.jday0;
TLm.iiAN_corr=iiAN_corr;
if iiAN_corr, TLm.AN_freqs=AN_freqs; end
%
% Add depth scale if desired:
if ii_geom,
    TLm.scales2=TLm.scales; TLm.labels2=TLm.labels;
    TLm.scales2{2}=xyz_all(iiph,3); TLm.labels2{2}='Depth - m';    
end
%
if ~isempty(out_root),
    if strcmp(suf,'_SL'),
        str=['SL=TLm;']; eval(str);
    	str=['save ',out_file,' SL']; eval(str)
    elseif strcmp(suf,'_AN'),
        str=['AN=TLm;']; eval(str);
    	str=['save ',out_file,' AN']; eval(str)
    else
    	str=['save ',out_file,' TLm']; eval(str)
    end
end
%
% COPYRIGHT (c) 2008
% 
% Applied Research Laboratories, The University of Texas at Austin
% (ARL:UT)
% All rights reserved.
% 
% This material may be reproduced by or for the U.S. Government
% only. This notice must appear in all copies of this file and its
% derivatives. Further distribution outside the U.S. Government is
% prohibited without the explicit written authorization from the
% Director of the Environmental Sciences Laboratory, ARL:UT.
% 
% THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
% EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE
% IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
% PURPOSE, OR NONINFRINGEMENT.
% 
% IN NO EVENT WILL ARL:UT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING
% OUT OF THE USE OR INABILITY TO USE THIS SOFTWARE, EVEN IF ARL:UT HAS
% BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
