mbs=bsModel(sgl, ts, cts, mch)

function mbs=bsModel(sgl, ts, cts, mch) 

%%----------------------------------------------------
%% Model slow baseline for TOF spectrum as exponential, Gaussian or linear
%%
%% INPUT:
%% sgl - observed unprocessed discrete TOF signal
%% ts - time samples (can be with unequal steps) for signal-less (noise) regions over
%%      the "decaying" part for baseline modeling 
%%      if "ts" set to a single number (length < 4), will call  "findBsEnv" for
%%      automated basline envelope estimate for modeling
%% cts - time interval for the estimate of constant offset
%% mch - choice of model: 1 - linear, 2 - exponential, 3 - Gaussian
%%
%% OUTPUT:
%% mbs - model baseline 
%%
%% NOTES:
%% --1-- Created by DM ([e-mail]) 10/10/04
%% --2-- No self-optimization code is included yet, but defaults work
%
% USAGE:
% mbs=bsModel(sgl, ts, cts, mch);
% Dependency: uses external "polyfit" function 
%

cfac = 0.05; %% correction for STD/mean of noise
rstp = 199; %% default down-sampling step for envelope finding;

sis=size(sgl); sit=size(ts); sict=size(cts);
if sis(1) < sis(2) sgl=sgl'; end;
if sit(1) < sit(2) ts=ts'; end;
if sict(1) < sict(2) cts=cts'; end;

cbs = mean(sgl(cts)); %% constant baseline offset
t1=1:length(sgl); t1=t1';
mbs = 0*sgl;

if length(ts) < 4  %% find baseline envelope points
    [ts, sts]=findBsEnv(sgl-cbs, rstp);
    sis1 = size(sts); sit1 = size(ts);
    if sis1(1) < sis1(2) sts=sts'; end;
    if sit1(1) < sit1(2) ts=ts'; end;
else %% use baseline envelope points provided by user
    sts = sgl(ts)-cbs; 
end;
cbs;
if (mch == 1) %% linear model
    [p1,s1]= polyfit(ts, sts, 1);
    mbs = p1(1)*t1+p1(2); 
    clear s1
elseif (mch == 2) %% exponential model
    sts(find(sts < 1e-3))=1; % prevent complex "log" 07/01/09
    [pe1,se1] = polyfit(ts, log(sts),1);
    mbs = exp(pe1(1)*t1+pe1(2));
    clear se1
else %% Gaussian model
    sts(find(sts < 1e-3))=1; % prevent complex "log" 07/01/09
    [pe2,se2] = polyfit(ts, log(sts),2);
    mbs = exp(pe2(1)*t1.*t1+pe2(2)*t1+pe2(3));
    [ymin,imin]=min(mbs);
    ii1=find(mbs(imin+1:end)>ymin);
    if length(ii1)>0 mbs(imin+ii1)=ymin; end;
    clear se2    
end;

mbs = (1-cfac)*mbs + cbs; %% correct by STD of noise to minimize negative signal

%'bsModel_debug:', sum(imag(mbs)), imag(cbs) %debug

return;