function X = simulateTimeresolved( t, Amplitudes, mySpectra, Tau, dBnoise, varargin)
% Simulate time resolved emission fluorescence (TRES) measurements. 
% Simulates the IRF (or uses measured one), decays with different lifetimes
% at different emission wavelengths, mixtures of these and a convolved
% version with noise (Data). 
%
%INPUT:
% t             A vector with time axis (in ns)
% Amplitudes    A matrix with amplitudes of the wanted fluorophore decays.
%               Each row should correspond to a "sample", so for example
%               Amplitudes=[1000 500 800; 100 600 1000] will have two
%               samples with a mixture of three fluorophores with different
%               amplitudes.
% mySpectra     A matrix with the spectral profiles of the fluorophores
%                Emission x Fluorophore 
% Tau           A vector of the fluophore lifetimes (in ns)
% dBnoise       Added homoscedastic noise with the poisson distribution
%
% Optional  
%               'measuredIRF': a vector with a measured IRF to use in stead
%                of the simulations.
%               'IRFparameters': A vector with the parameters 
%                [cIRF, k1IRF,k2IRF] for the instrument response function
%                (IRF). Default is [1000, 15, 14.5]
%
%OUTPUT:
% X             Struct with the simulated Data (time decay x samples x
%               emissions), the IRF, the timeaxis (t), etc.
%               The Pure, Mixed and Convoluted data is in the "Extra" field. 

% Based on the simulations described by: https://link.springer.com/article/10.1007/s10895-009-0589-1
% Copyright, 2023 - 
% This M-file and the code in it belongs to the holder of the
% copyrights and is made public under the following constraints:
% It must not be changed or modified and code cannot be added.
% The file must be regarded as read-only. 
%
% ?smund Rinnan
% Department of Food Science
% Faculty of Sciences
% University of Copenhagen
% Rolighedsvej 26, 1958 Frederiksberg C, Denmark


% Optional inputs
i = 1;
if (~isempty(varargin))
    while i <= size(varargin,2)
        switch varargin{i}
            case {'measuredIRF'}
                measuredIRF= varargin{ i + 1};
                X.IRF = measuredIRF;
                i = i + 1;
            case {'IRFparameters'}
                IRFparameters= varargin{i+1};
                if length(IRFparameters) == 3
                    cIRF = IRFparameters(1);
                    k1IRF = IRFparameters(2);
                    k2IRF = IRFparameters(3);
                else
                    error('Simulated IRF must be a vector with 3 elements')
                end
                i = i + 1;
            otherwise
                error(['Invalid optional argument'])
        end
        i = i + 1;
    end
elseif isempty(varargin) % check if optional inputs about IRF have been set, otherwise set to default. 
    cIRF = 10000000;
    k1IRF = 15;
    k2IRF = 14.5;
end

if size( mySpectra, 2) ~= length( Tau)
    error( 'The number of fluorophores in the pure spectra, and Tau should be the same!')
end

X.t=t;
X.details.Lambda = dBnoise;
X.Extra.PureSpectra = mySpectra;

vec = @(t) t(:);

%% simulate instrument response function: 

if ~isfield(X,'IRF') % only do it if we are not using a measured one. 
    X.IRF = cIRF*t.^2.*((exp(-k2IRF*t)-exp(-k1IRF*t)));
    X.details.cIRF = cIRF;
    X.details.k1IRF = k1IRF;
    X.details.k2IRF = k2IRF;
end

%% simulate chromphore decays:

% to collect the mixed decay
X.Extra.Mixed=zeros( length(t), size(Amplitudes,1), size(mySpectra,1)); 

X.Amplitudes = Amplitudes;
X.Tau = Tau;
X.Extra.PureExp = (exp(-t./Tau));

for n=1:size( Amplitudes, 1)
    for j=1:size( mySpectra, 1)
        X.Extra.PureDecay( :, :, n, j) = Amplitudes( n, :) .* mySpectra( j, :) .* exp(-t./Tau);
        X.Extra.Mixed( :, n, j) = X.Extra.Mixed( :, n, j) + sum( X.Extra.PureDecay( :, :, n, j),2);
    end
end

%% Convolve with the IRF:

X.IRF = normalize( X.IRF(:));
for n=1:size(Amplitudes,1)
    for j=1:size(mySpectra,1)
        Convoluted = conv( X.Extra.Mixed( :, n, j), X.IRF);
        X.Extra.Convoluted( :, n, j) = round( Convoluted( 1:length(t)));
    end
end

%% Add noise

vecsig = vec( X.Extra.Convoluted);
sig = sqrt( mean( vecsig.^2));

mynoise = poissrnd( 1e5, size( X.Extra.Convoluted));
%Set the smallest noise to 0 (only want positive noise)
mynoise = mynoise - min( mynoise(:));
vecnoise = vec( mynoise);
%mynoise needs to be scaled
mynoise = mynoise/ sqrt( mean( vecnoise.^2));
finalnoise = mynoise * sig * 10^(-dBnoise/10);
X.Data = X.Extra.Convoluted + finalnoise;

