function TSP=orca_tspread(SD,f_lims,zs_lims,zr_lims,rng_lims,tau_lims,BW,iinorm,tsize,iifig,iinewfig);
% Call TSP=orca_tspread(SD,f_lims,zs_lims,zr_lims,rng_lims,tau_lims,BW,iinorm,tsize,iifig,iinewfig);
%
%
% 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.
if nargin<8, iinorm=1; end
if nargin<9, tsize=[]; end
if nargin<10, iifig=1; end
if nargin<11, iinewfig=1; end
%
% Read in needed
jset_phi=find(strcmp(SD.data_name,'|phi|'));
if isempty(jset_phi),
    disp('|phi| not found');
end
%
[PHIR,stat]=hdfsd_read(SD,{zr_lims,[],f_lims},jset_phi); zr=PHIR.scales{1}; nzr=length(zr);
[PHIS,stat]=hdfsd_read(SD,{zs_lims,[],f_lims},jset_phi); zs=PHIS.scales{1}; nzs=length(zs);
freq=PHIR.scales{3}; nfreq=length(freq);
%
jset_knr=find(strcmp(SD.data_name,'Re(kn)'));
[KNR,stat]=hdfsd_read(SD,{f_lims,[]},jset_knr);
jset_kni=find(strcmp(SD.data_name,'Im(kn)'));
[KNI,stat]=hdfsd_read(SD,{f_lims,[]},jset_kni);
jset_vg=find(strcmp(SD.data_name,'Group Speed - m/s'));
[VG,stat]=hdfsd_read(SD,{f_lims,[]},jset_vg);
%
rng_km=adi_lim_eval(rng_lims); nrng=length(rng_km); rng=1000*rng_km;
%
Ts_ax=spread_env_init2(tau_lims,1/BW,3);
tau=Ts_ax.xax;
ntau=Ts_ax.nx;
fac=sqrt(2*pi);
TSP.file_name=SD.hdf_file;
TSP.dlab='TL - dB';
TSP.data=zeros(nzr,nrng,ntau,nzs,nfreq);
TSP.ndim=5;
TSP.dims=[nzr,nrng,ntau,nzs,nfreq];
TSP.labels={'Depth - m','Range - km','Time Delay - s','Source Depth - m','Frequency - Hz'};
TSP.scales={zr,rng_km,tau,zs,freq};
TSP.iidB=1;
if iinorm,
    TSP.dB_max=zeros(nrng,nfreq);
end
TSP.vg_ref=zeros(1,nfreq);
%
for jf=1:nfreq,
    iim=find(~isnan(VG.data(jf,:))); nmode=length(iim);
    [im_min,jm_ref]=min(KNI.data(jf,iim)); vg_ref=VG.data(jf,jm_ref);
    tau_norm=rng/vg_ref; TSP.vg_ref=vg_ref;
	for jzs=1:nzs,
		for jzr=1:nzr,
            phi_fac=fac*PHIS.data(jzs,iim,jf).*PHIR.data(jzr,iim,jf);
			for jrng=1:nrng,
				taux=rng(jrng)./VG.data(jf,iim) - tau_norm(jrng);
				p_mode=phi_fac.*exp(-KNI.data(jf,iim)*rng(jrng))./sqrt(KNR.data(jf,iim)*rng(jrng));
				psq_mode=max(1e-30,p_mode.^2);
				TSP.data(jzr,jrng,:,jzs,jf)=10*log10(spread_env2(Ts_ax,taux,psq_mode)');
			end
		end
	end
    if iinorm,
        for jrng=1:nrng,
            dB_max=max(max(max(TSP.data(:,jrng,:,jzs,jf),[],3),[],1),[],4);
            TSP.data(:,jrng,:,:,jf)=TSP.data(:,jrng,:,:,jf) - repmat(dB_max,[nzr,1,ntau,nzs]);
            TSP.dB_max(jrng,jf)=dB_max;
        end
    end
end
dBmax=ceil(max(TSP.data(:))/5)*5;
TSP.clim=[dBmax-40,dBmax];
%
if iifig,
    [dim_sort,ii]=sort([nzr,nrng,0,nzs,nfreq]); vdims=[ii(5),3,ii(4)],
    view5(TSP,vdims,[],iinewfig,tsize);
end
