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

%Instrument response function convolution with data (according to BBA van Stokkum 2004
%1657 82-104)
%Dispersion function could be included, as in van Global and target analysis of time-
%resolved spectra Lecture notes Troisième Cycle de la Physique en Suisse Romande, March
%14-24, 2005 Ivo H.M. van Stokkum
%This function is called by globfun_bio
%
%IRF_fun is the generated IRF concentration profile

function [decays,IRF_fun] = glob_fun_IRF(start_pars,start_pars_nonlin,Par_fix_time_mat,decays,time_fix_counter)

global Par_fix_IRF comps Par_fix_time_zero xdata rev_counter

fix_counter=0; %Count fixed time constants again as IRF is built up
if Par_fix_IRF==0 %If IRF parameters need to be fitted
    FWHM_gauss_tilde=start_pars_nonlin(comps+2-time_fix_counter-Par_fix_time_zero)/(2*sqrt(2*log(2)));
    for i=1:comps-rev_counter
        if isequal(i,Par_fix_time_mat(i))==0
            IRF_comps(:,i)=0.5*exp(start_pars_nonlin(i-fix_counter)*(start_pars_nonlin(comps+3-time_fix_counter-Par_fix_time_zero)+0.5*start_pars_nonlin(i-fix_counter)*FWHM_gauss_tilde^2))*...
                (1+erf((xdata-(start_pars_nonlin(comps+3-time_fix_counter-Par_fix_time_zero)+start_pars_nonlin(i-fix_counter)*FWHM_gauss_tilde^2))/sqrt(2)/FWHM_gauss_tilde));
        elseif isequal(i,Par_fix_time_mat(i))==1
            IRF_comps(:,i)=0.5*exp(start_pars(i)*(start_pars_nonlin(comps+3-time_fix_counter-Par_fix_time_zero)+0.5*start_pars(i)*FWHM_gauss_tilde^2))*...
                (1+erf((xdata-(start_pars_nonlin(comps+3-time_fix_counter-Par_fix_time_zero)+start_pars(i)*FWHM_gauss_tilde^2))/sqrt(2)/FWHM_gauss_tilde));
            fix_counter=fix_counter+1;
        end
        decays_IRF(:,i)=decays(:,i).*IRF_comps(:,i);
    end
    clear i fix_counter
    IRF_fun=(1/(FWHM_gauss_tilde*sqrt(2*pi)))*exp(-log(2)*(2*(xdata-start_pars_nonlin(comps+3-time_fix_counter-Par_fix_time_zero))/start_pars_nonlin(comps+2-time_fix_counter-Par_fix_time_zero)).^2);
    decays=decays_IRF;
elseif Par_fix_IRF==1 %If IRF parameters are fixed
    FWHM_gauss_tilde=start_pars(comps+2)/(2*sqrt(2*log(2)));
    for i=1:comps-rev_counter
        if isequal(i,Par_fix_time_mat(i))==0
            IRF_comps(:,i)=0.5*exp(start_pars_nonlin(i-fix_counter)*(start_pars(comps+3)+0.5*start_pars_nonlin(i-fix_counter)*FWHM_gauss_tilde^2))*...
                (1+erf((xdata-(start_pars(comps+3)+start_pars_nonlin(i-fix_counter)*FWHM_gauss_tilde^2))/sqrt(2)/FWHM_gauss_tilde));
        elseif isequal(i,Par_fix_time_mat(i))==1
            IRF_comps(:,i)=0.5*exp(start_pars(i)*(start_pars(comps+3)+0.5*start_pars(i)*FWHM_gauss_tilde^2))*...
                (1+erf((xdata-(start_pars(comps+3)+start_pars(i)*FWHM_gauss_tilde^2))/sqrt(2)/FWHM_gauss_tilde));
            fix_counter=fix_counter+1;
        end
        decays_IRF(:,i)=decays(:,i).*IRF_comps(:,i);
    end
    clear i fix_counter
    IRF_fun=(1/(FWHM_gauss_tilde*sqrt(2*pi)))*exp(-log(2)*(2*(xdata-start_pars(comps+3))/start_pars(comps+2)).^2);
    decays=decays_IRF;
end
clear decays_IRF