doRTReaction(tmpl,dn,pp,cs,varargin)

function [rn,ctpoint,ct] = doRTReaction(tmpl,dn,pp,cs,varargin)
%DOTRIPLERT(tmpl,dn,pp,cs,[eff,makeplot,nf,bf,col]);
%
% This function will create amplification (Rn) data from an artificial
% amount of cDNA and it will calculate the Ct value of the reaction.  The
% treshold to calculate the Ct value will be placed on Rn = 1.
%
%INPUT PARAMETERS IN THE CORRECT ORDER, SQUARE BRACKETS MEAN THAT THE
%PARAMETER IS OPTIONAL
% - Template (tmpl): amount of added cDNA (total RNA!) template in nanograms
% - DNTP's (dn): amount of added DNTP's also in nanograms
% - Primers & Probes (pp): amount of primers and probes
% - Cycles (cs): number of replication cycles
% - [Efficiency (eff)]: efficiency of the primer-probe set (0-2), if omitted E = 1.8
% - [Makeplot(mp)]: make an amplification plot, not if omitted
% - [Noising factor (nf)]: set an artificial noising factor, if omitted 
%                          NF = 0.05, this noising factor is actually the   
%                          machine related detection error of the 
%                          fluorescent light intensity
% - [Bleaching factor (bf)]: set the fluorescence bleaching factor, if omitted
%                            BF = 0.0001. Every cycle, the fluorescent reaction
%                            products are exited, inducing bleaching after a
%                            number of cycles.  This is a minor effect that 
%                            becomes important after a large number of cycles.
%                            Stable fluorophores have a low bleaching factor,
%                            unstable ones a high factor, eg 0.005.
% - [Color (col)]: Color of the polotted line (if Makeplot = 1), by 
%                  default blue.  If no color is specified, a random 
%                  color is generated.
%
%OUTPUT PARAMETERS
% - Rn: returns the datapoints of the amplification curve
% - CtPoint: the cordinates of the Ct point
% - Ct: Ct value
%
% (c)2007 Gie Spaepen
% University of Antwerp faculty of Medicine
% Universiteitsplein 1, Building T3.33
% 2610 Wilrijk
%
% Comments: [e-mail]

% Program vars
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Flag if a plot has to be drawn
makeplot = 0;
% Default efficiency
eff = 1.8;
% Noising factor; see help above
nf = 0.05;
% Fluorescence bleaching factor, see help above
bf = 0.001;
% Default line color
col = 'b';
% Resultmatrix
results = [1:cs];
% Cyclematrix for plotting
csm = [1:cs];

% Adjust parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The initial amount is given in nanograms, but
% this is total RNA, but the real number of mRNA
% copies is much lower.  To generate a realistic
% PCR effect we devide the amount of template by
% 10e6.
tmpl = tmpl/1000000;

% Check input parameters...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ...on efficiency
if nargin >= 5
    eff = varargin{1};
end
% ...if a plot should be drawn
if nargin >= 6
    makeplot = varargin{2};
end
% ...if another noising factor should be used
if nargin >= 7
    nf = varargin{3};
end
% ...if another bleaching factor should be used
if nargin >= 8
    bf = varargin{4};
end
% ...if you choose a color for the plot
if nargin >= 9
    col = varargin{5};
end
    
% Intermediate/Temp vars
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ... for the amount of DNTP's
t_dn = dn;
% ... for the amount of primers/probes
t_pp = pp;
% ... for the amount of cDNA
t_tmpl = tmpl;
% ... for the fluorescent signal after every cycle
t_rn = 0;
% ... the bleaching factor after every cycle
t_bf = 0;

% Do reaction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:cs
    % Calculate the % of used DNTP's
    t_dn_eff = t_dn/dn;
    % Calculate the % of used primers/probes...
    t_pp_eff = t_pp/pp;
    % Calculate the efficiency
    t_eff = eff * t_dn_eff * t_pp_eff;
    % Get the reaction product
    t_p = t_tmpl*(t_eff^2);
    % Update the temporal variables after reaction by substracting the 
    % amount of reaction product from the pool of primers and DNTP's
    t_tmpl = t_p;
    t_dn = t_dn - t_p;
    t_pp = t_pp - t_p;
    % Add the fluorescence plus noising factor
    t_rn = t_rn+(t_p+((rand(1)*nf)-(nf/2)));
    % Substract the fluorescence bleaching (minor effect)
    t_bf = (bf*((2*i)+1));
    % If the bleaching effect is smaller than the fluorescence 
    % the bleaching is substracted.  Otherwise, we are close to 
    % the detection point (0) and nothing will be done.
    if t_bf < t_rn
        t_rn = t_rn - t_bf;
    end  
    % Update the matrix
    results(i) = t_rn;
end

% Calculate Ct Value by finding the 2 points, one just above and one just
% below the treshold
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
try
    % Get the X values from the result matrix
    x2 = find(results > 1,1);
    x1 = x2-1;
    % Get the Y values from the result matrix
    y1 = results(x1);
    y2 = results(x2);
    % Calculate the line parameters (a,b) of these 2 points 
    ab = [x1,1;x2,1]\[y1;y2];
    a = ab(1);
    b = ab(2);
    % The treshold is set to 1 so this line is treshold -> Y = 1.  The crossing
    % of the two lines is the Ct point.  ct is then ctpoint(2)
    ctpoint = [1,-a;1,0]\[b;1];
    ct = ctpoint(2);
catch
    % The values cannot be detected so we return a 0 ct
    disp('* Ct cannot be detected 0 values returned');
    ctpoint = [0 0];
    ct = 0;
end


% Make plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if makeplot == 1
    % Plot it
    plot(csm,results,col);
    % Set the labels
    xlabel('Cycles');
    ylabel('dRn');
    hold on;
    % Plot the points
    plot(ctpoint(2),ctpoint(1),strcat(col,'o'));
    hold off;
end

% Print results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(strcat('* Ct  : ',num2str(ct)));

rn = results;