function [T_peaks P_peaks] = Find_peaks_TP(signal,Q_peaks,R_peaks,S_peaks,I,status)
% [T_peaks P_peaks] =
% Find_peaks_TP(signal,Q_peaks,R_peaks,S_peaks,I,status)
%
% Input Variables
%   signal  - The difference between two cell voltages
%   Q_peaks - The set of Q peaks
%   R_peaks - The set of R peaks
%   S_peaks - The set of S peaks
%   I       - The set of peaks        
%   status  - the status of peaks (min or max)
% Output Variables
%   T_peaks - The set of T peaks
%   P_peaks - The set of P peaks
%
% Description
%   Determines the T and P wave peaks
% Created
%   8/23/2006, Mohamed A Mneimneh, Marquette University
%
% This software is released under the terms of the GNU General
% Public License (http://www.gnu.org/copyleft/gpl.html).
% Modified

count = 1;
T_peaks=1;
P_peaks=1;
for i = 1:length(R_peaks)-1
    X1 = 0;
    X2 = 0;
    x_p = 0;
    flag_p = 0;
    flag_t =0;
    x_t =0;
    k1 = S_peaks(i) + floor(0.14*(R_peaks(i+1)-R_peaks(i)));
    k2 = Q_peaks(i+1)-10;
    clear intercept
    for j = floor(abs(signal(R_peaks(1)))):-10:0

        intercept = j*ones(k2-k1+1,1);
         if ( length(X1) < 4 )
            X1 = k1 + intersect(intercept, signal(k1:k2));
            if length(X1) == 3
                if ( (X1(2)- X1(1)) > (X1(3)- X1(2)) )
                    X1 = X1(2:3);
                    k1 = X1(1) +10;
                    intercept = j*ones(k2-k1+1,1);
                    X1 = k1 + intersect(intercept, signal(k1:k2));
                else
                    k2 = X1(3) -10;
                    X1 = X1(1:2);
                    intercept = j*ones(k2-k1+1,1);
                    X1 = k1 + intersect(intercept, signal(k1:k2));
                end
            end
            if (length(X1) > 1 )
                if ( ( abs(X1(1) - S_peaks(i)) < abs(X1(2) -  Q_peaks(i+1)) || flag_p) && ~flag_t)
                    k1 = X1(2)+ floor(0.4*abs(X1(2) -k2));
                    intercept = j*ones(k2-k1+1,1);
                    x_t = X1;
                    flag_t = 1;
                elseif ~flag_p && length((status(find( I == [find(signal == max(abs(signal(X1(1) :X1(2))))) find(signal == -max(abs(signal(X1(1) :X1(2)))))])))) && (status(find( I == [find(signal == max(abs(signal(X1(1) :X1(2))))) find(signal == -max(abs(signal(X1(1) :X1(2)))))])))
                    k2 = S_peaks(i)+ floor(0.8*abs(X1(1) -k1));
                    intercept = j*ones(k2-k1+1,1);
                    x_p = X1;
                    flag_p = 1;
                end
            end
        end
        if ( length(X2) < 4 )
            X2 =k1 + intersect(-intercept, signal(k1:k2));
            if length(X2) == 3
                if ( (X2(2)- X2(1)) > (X2(3)- X2(2)) )
                    X2 = X2(2:3);
                    k1 = X2(1) +10;
                    intercept = j*ones(k2-k1+1,1);
                    X2 =k1 + intersect(-intercept, signal(k1:k2));
                else
                    k2 = X2(3)-10;
                    X2 = X2(1:2);
                    intercept = j*ones(k2-k1+1,1);
                    X2 =k1 + intersect(-intercept, signal(k1:k2));
                end
            end
            if ( length(X2) > 1 ) 
                if ( ( abs(X2(1) - S_peaks(i)) < abs(X2(2) -  Q_peaks(i+1)) || flag_p ) && length(x_t) == 1 && ~flag_t)
                    k1 = X2(2)+ floor(0.4*abs(X2(2) -k2));
                    intercept = j*ones(k2-k1+1,1);
                    x_t = X2;
                    flag_t = 1;
                    
                elseif ~flag_p && length((status(find( I == [find(signal == max(abs(signal(X2(1) :X2(2))))) find(signal == -max(abs(signal(X2(1) :X2(2)))))])))) > 1 && (status(find( I == [find(signal == max(abs(signal(X2(1) :X2(2))))) find(signal == -max(abs(signal(X2(1) :X2(2)))))])))
                
                    k2 = S_peaks(i)+ floor(0.8*abs(X2(1) -k1));
                    intercept = j*ones(k2-k1+1,1);
                    x_p = X2;
                    flag_p = 1;
                end
            end
            X2 =0;
        end


        if ( length(x_p) == 2 )
            P_peaks(count) = [find(signal == max(abs(signal(x_p(1) :x_p(2))))) find(signal == -max(abs(signal(x_p(1) :x_p(2)))))];
          end
        if (length(x_t) == 2)
            T_peaks(count) = [find(signal == max(abs(signal(x_t(1) :x_t(2))))) find(signal == -max(abs(signal(x_t(1) :x_t(2)))))];
        end
        if (length(x_t) == 2 && length(x_p) == 2)
            count = count +1;
            clear X1 X2 x_p x_t flag_p flag_t
            break;
        end
    end
end