%Copyright (c) 2010, Luuk J.G.W. van Wilderen
%
%globfun_spectra.m

function [spectra]=globfun_spectra(decays,start_pars_lin,rev_counter)
%Calculation of spectra.
%Fitting of linear parameters by Matlab built-in function of mlr (matrix
%left division, mldivide), which is very fast.
%If the errors need to be calculated, lscov is used. Both mldivide and lscov
%use QR decomposition, and give identical results. However 'lscov' is
%slightly slower, so therefore this is only used when the errors are
%requested. It is also possible that mldivide convergences to a slightly
%less accurate solution, but these difference are minimal (of the order of 1e-4).
%NB. Last decay curve not fitted, as its spectrum is zero! This is because
%the before-last species needs to decay into another species (without
%amplitude).

global Par_IRF Par_fix_OD_offset Error_switch_lin_spectra 
global comps opt_pars_lin opt_pars_lin_std spectra_std ydata
% %Defined in dispersion_IRF_calc
global delay_shift_curve xdata dispersion_correction

%Dispersion correction applied to data (ydata). This script does not modify
%the data, but only corrects if for correct spectrum calculation.
%For every channel, each time point is shifted with a specific delay time
%as determined by the dispersion curve.
if dispersion_correction==1
    for i=1:size(ydata,2) %Fast interpolator chosen because time points are most likely irregularly spaced
        yi = interp1q(xdata,ydata(:,i),xdata+delay_shift_curve(i));
        ydata_disp(:,i)=yi;
    end
    ydata_disp(isnan(ydata_disp)) = 0;  %Remove any NaN's ('interp1q' above cannot extrapolate, so leaves NaN's)
elseif dispersion_correction==0         %The fitting function can't handle NaN's
    ydata_disp=ydata;
end

if Error_switch_lin_spectra == 0
    if Par_fix_OD_offset == 0
        spectra_temp=(decays(:,1:comps+Par_IRF-rev_counter))\ydata_disp; 
        res_temp=sum(sum(decays(:,1:comps+Par_IRF-rev_counter)*spectra_temp-ydata_disp).^2);         %OD-offset is determined by fitting a zeroth order polynomial to residual
        [opt_pars_lin,opt_pars_lin_err]=polyfit(ones(size(res_temp)),res_temp,0);   %Linear fit for baseline (zeroth order)
        opt_pars_lin_std=opt_pars_lin_err.normr;
        y_fit_temp=ydata_disp+opt_pars_lin;
        spectra=decays(:,1:comps+Par_IRF-rev_counter)\y_fit_temp;                               %Added baseline. Last decay curve not fitted, as its spectrum is zero (Last species decays in species with zero amplitude)
        clear spectra_temp y_fit_temp
    elseif Par_fix_OD_offset == 1
        opt_pars_lin=start_pars_lin;
        y_fit_temp=ydata_disp+opt_pars_lin;
        spectra=decays(:,1:comps+Par_IRF-rev_counter)\y_fit_temp;   
        clear y_fit_temp
    end
elseif Error_switch_lin_spectra == 1 %Because errors are requested, lscov is used instead of '\' (mldivide). Is slightly slower than '\', so that is default.
    if Par_fix_OD_offset == 0
        [spectra_temp]=lscov((decays(:,1:comps+Par_IRF-rev_counter)),ydata_disp);
        res_temp=sum(sum(decays(:,1:comps+Par_IRF-rev_counter)*spectra_temp-ydata_disp).^2);         %OD-offset is determined by fitting a zeroth order polynomial to residual
        [opt_pars_lin,opt_pars_lin_err]=polyfit(ones(size(res_temp)),res_temp,0);   %Linear fit for baseline (zeroth order)
        opt_pars_lin_std=opt_pars_lin_err.normr;
        y_fit_temp=ydata_disp+opt_pars_lin;
        [spectra,spectra_std]=lscov((decays(:,1:comps+Par_IRF-rev_counter)),y_fit_temp);
        clear spectra_temp y_fit_temp
    elseif Par_fix_OD_offset == 1
        opt_pars_lin=start_pars_lin;
        y_fit_temp=ydata_disp+opt_pars_lin;
        [spectra,spectra_std]=lscov((decays(:,1:comps+Par_IRF-rev_counter)),y_fit_temp); 
        clear y_fit_temp
    end
end