function [CTS,rspec,FFT2]=cts2spec2(cts_root,FFT,fax_get,fzpad,T_sec,Df)
%
% Differs from cts2spec in that a frequency shift Df may be applied to the replica first;
%
% Computes a spectrum rspec from a source waveform in CTS file with cts_root that matches
% the spectra in the FFT file given by the structure FFT (from fft_read1_open).  The 
% frequency band to be processed is given by fax_get.
%
% See:  cts_details, cts_open, cts_read, cts_write, cts_write_open, cts_combine, cts2fft, cts2spec, 
%       cts_from_real, cts_shot, cts_wf_make, fft2cts, mcs_gui.
%
%
% 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<5, T_sec=0; end
if nargin<6, Df=0; end
%
CTS=cts_open(cts_root); f0=CTS.f_bb_shift;
f_lims=[fax_get(1), fax_get(end)]; nfreq=length(fax_get);
if prod(f_lims-f0)>0,
    disp(f0); disp(f_lims);
    error('CTS file center frequency not inside f_lims: ');
end
dop_alpha=(Df + f0)/f0;
sig_bb=cts_read(CTS,0,inf); npt=length(sig_bb); FFT2.PRI=npt/CTS.fs/dop_alpha; df_bb=1/FFT2.PRI;
%
% Match df so that time window same:
FFT2.df=FFT.df;  FFT2.Tw=1/FFT.df;
% Make FFT2.fs such that FFT2.dt fits by power-of-two integer into FFT2.Tw and FFT2.fs>fband:
FFT2.nfft=2^nextpow2(diff(f_lims)/FFT2.df);

if fzpad<0, disp('fzpad should now be >=0 for enveloping also'); fzpad=1; end
fband=diff(f_lims);
fac=2^(-(nextpow2(FFT.fs/fband)-1-fzpad)); fsx=FFT.fs*fac; disp(['fsx from cts2spec = ',num2str(fsx)]);
%
FFT2.nfft=FFT.nfft*fac; FFT2.nT_int=fac*FFT.nT_int;
FFT2.dt=FFT2.Tw/FFT2.nfft; FFT2.fs=1/FFT2.dt;
dt_new=dop_alpha*FFT2.dt;
npt_rep=round(FFT2.PRI/dt_new), old_npt=round(FFT2.PRI/FFT2.dt),
%
t_bb=[0:npt]'/CTS.fs; 
s_bb=[sig_bb(:); sig_bb(1)];
t_uni=[0:npt_rep-1]' * FFT2.dt; 
p_rep0=zeros(FFT2.nfft,1); 
p_rep0(1:npt_rep)=interp1(t_bb,s_bb,dop_alpha*t_uni).*exp(i*2*pi*Df*t_uni);
%
FFT2.f_bb_shift=round(f0/FFT2.df)*FFT.df;
%delf=f0 - FFT2.f_bb_shift;     % DO NOT CHANGE f_bb_shift of CTS file!!!
%jf_get1=round((fax_get(1)-FFT2.f_bb_shift)/FFT.df) + 1;
%if jf_get1<=0, jf_get1=FFT2.nfft+jf_get1; end
%iineg=[jf_get1:FFT2.nfft]; lneg=length(iineg);
%jf_get2=round((fax_get(end)-FFT2.f_bb_shift)/FFT.df) + 1;
%iipos=[1:jf_get2];

jf_bb=round((FFT2.f_bb_shift - fax_get(1))/FFT2.df) + 1;
iipos2=[jf_bb:nfreq]; npos2=length(iipos2); iiS_pos=[1:npos2];
iineg2=[1:jf_bb-1]; iiS_neg=FFT2.nfft + [-jf_bb+2:0];
%
% Compute number of T_sec sections of waveform to split signal into:
if T_sec==-999,
    FFT2.ndt_sec=FFT2.nfft;
	FFT2.n_sec=1;
	rspec=zeros(nfreq,1);
    npri=ceil(FFT2.Tw/FFT2.PRI),
    p_rep=repmat(p_rep0(1:npt_rep),[1,npri]);
    p_rep=p_rep(1:FFT2.nfft);
    S_bb=fft(p_rep);
%	rspec(1:lneg)=S_bb(iineg);
%	rspec(lneg+1:nfreq)=S_bb(iipos);
    rspec(iipos2)=S_bb(iiS_pos);
    rspec(iineg2)=S_bb(iiS_neg);
else
    if T_sec>0, n_sec=round(FFT2.PRI/T_sec); else, n_sec=1; end
	FFT2.ndt_sec=round(npt_rep/n_sec);
	FFT2.n_sec=n_sec;
%
	rspec=zeros(nfreq,n_sec);
%
	for j_sec=1:n_sec,
        p_rep=zeros(size(p_rep0));
        ii=(j_sec-1)*FFT2.ndt_sec + [1:FFT2.ndt_sec];     % Non-zero portion of section j_sec;
% This puts section in its usual place in the original waveform:
%        p_rep(ii)=p_rep0(ii);
% This puts each section at the beginning of the waveform:
        p_rep(1:FFT2.ndt_sec)=p_rep0(ii);
        S_bb=fft(p_rep);
%		rspec(1:lneg,j_sec)=S_bb(iineg);
%		rspec(lneg+1:nfreq,j_sec)=S_bb(iipos);
        rspec(iipos2,j_sec)=S_bb(iiS_pos);
        rspec(iineg2,j_sec)=S_bb(iiS_neg);
	end
end
