function [QRS1,Rp,Sp,R2p,QRS2,ymaxaux,type,Sgran]=Rwave(n,X,Xpb,D,Der,PKni,M,I,Fs,Kr,Ks,Krr)

% ---- QRS complex classification ----
% In RSR' type significative points are returned,
% in QRS and QR types, R and S wave positions, are returned respectively.

% ---- Initialization ----
QRS1=[]; QRS2=[]; Qp=[]; Rp=[]; Sp=[]; R2p=[]; Rex=0; Qex=0; Sex=0; R2ex=0;
type=0; noR=0; Sgran=0;

% ---- Previous and former peaks ----

Ir=find(I>PKni); 
mpicd=I(Ir(1));
Ir=find(I<PKni);
mpici=I(Ir(length(Ir)));

% ---- RSR' type? ----
ydi=D(mpici);
ydd=D(mpicd);
ymaxaux=max(abs(ydd),abs(ydi));
kpi=2;

if (Xpb(PKni)<0)|(ydi>0&ydd<0&(kpi*ydi<(-1)*ydd|kpi*(-1)*ydd<ydi))
    perc=0.25;
    if (Xpb(PKni)>0&ydi>0&ydd<0)|((1+perc)*(-1)*ydi>ydd&(1-perc)*(-1)*ydi<ydd)
        % ---- RSR'type ----
        type=2;
        perc=0.35;
        if (Xpb(PKni)<0)
            % ---- PKni corresponds to S wave, R' will be on the right and R
            % on the left ----
            Daux=D(mpicd:length(D));
            ncero=zerocross(Daux);
            if isempty(ncero) return; end
            ncero=mpicd+ncero-1;
            Ir=find(I>ncero);
            if ~isempty(Ir) mpda=I(Ir(1));
                if ((-1*D(mpda)<ydd/5) & (abs(Xpb(ncero))<abs(Xpb(PKni))/10))  % JGM
                    type=3; % ---- Very large Q or S wave ----
                else
                    R2p=ncero;
                    Sp=PKni;
                    Daux=flipud(D(1:mpici));
                    ncero=zerocross(Daux);
                    if isempty(ncero) return; end
                    ncero=mpici-ncero+1;
                    Ir=find(I<ncero);
                    if ~isempty(Ir)
                        mpda=I(Ir(length(Ir)));
                        if ((D(mpda)<-ydi/5) & (abs(Xpb(ncero))<abs(Xpb(PKni))/10))  % JGM
                            type=3; % ---- Very large Q or S wave ----
                        else
                            Rp=ncero;
                        end
                    else
                         Rp=ncero;  % modified by guangyubin
                    end
                end
            else
                return ;
            end
        elseif abs(ydi)<abs(ydd)
            % ---- PKni corresponds to R wave ----
            Daux=D(mpicd:length(D));
            ncero=zerocross(Daux);
            if isempty(ncero) return; end
            ncero=mpicd+ncero-1;
            Ir=find(I>ncero);
            if ~isempty(Ir)
                mpic=I(Ir(1));
                if ~((1+perc)*abs(ydd)>abs(D(mpic))&(1-perc)*abs(ydd)<abs(D(mpic)))
                    type=1; % ---- Normal QRS type ----
                else
                    Sp=ncero;
                    Daux=D(Sp:length(D));
                    ncero=zerocross(Daux);
                    if isempty(ncero) return; end
                    ncero=Sp+ncero-1;
                    R2p=ncero;
                    Rp=PKni;
                end
            else
                return ;
            end
            
        elseif abs(ydi)>abs(ydd)
            % ---- PKnii corresponds to R' wave ----
            Daux=flipud(D(1:mpici));
            ncero=zerocross(Daux);
            ncero=mpici-ncero+1;
            if isempty(ncero)   % raw: if ~isempty(ncero)  modified by guangubin@bjut.edu.cn
                return;
            end
            Ir=find(I<ncero);
            if ~isempty(Ir)
                mpic=I(Ir(length(Ir)));
                if (~((1+perc)*abs(ydi)>abs(D(mpic))&(1-perc)*abs(ydi)<abs(D(mpic))))
                    type=1; % ---- Normal QRS type ----
                else
                    Sp=ncero;
                    Daux=flipud(D(1:Sp));
                    ncero=zerocross(Daux);
                    if isempty(ncero) return; end
                    ncero=Sp-ncero+1;
                    Rp=ncero;
                    R2p=PKni;
                end
            else
                return;
            end
        end
        if (type==2)&(R2p-Rp)/Fs>150e-3
            if Xpb(PKni)>0
                type=1;  % ---- Normal QRS type ----
            else type=3; % ---- Very large Q or S wave ----
            end
        end
    else type=3; % ---- Very large Q or S wave ----
    end
else type=1;   % ---- Normal QRS type ----
end


% ---- Onset and offset of RSR' ----
if (type==2)
    Ir=find(I>R2p);
     if isempty(Ir) return ;end;
     mpicd=I(Ir(1));
    Ir=find(I<Rp); 
    if isempty(Ir) return ;end;
    mpici=I(Ir(length(Ir)));
    R2p=testpeak(X,R2p,Fs,1);
    Sp=testpeak(X,Sp,Fs,0);
    Rp=testpeak(X,Rp,Fs,1);
    umbral=D(mpicd)/Krr;
    Daux=D(mpicd:length(D));
    QRS2=thresholdcross(Daux,umbral);
    QRS2=mpicd+QRS2-1;
    if isempty(QRS2) return; end
    umbral=X(QRS2);
    Xaux=flipud(X(1:R2p));
    S1=thresholdcross(Xaux,umbral);
    S1=R2p-S1+1;
    umbral=D(mpici)/Kr;
    Daux=flipud(D(1:mpici));
    QRS1=thresholdcross(Daux,umbral);
    QRS1=mpici-QRS1+1;
    if isempty(QRS1) return; end
    umbral=X(QRS1);
    Xaux=X(Rp:length(X));
    Q2=thresholdcross(Xaux,umbral);
    Q2=Rp+Q2-1;
end


% ---- R wave location in the very large Q or S wave case ----
if (type==3)
    R2p=[]; Rp=[]; Sp=[];
    Daux=D(mpicd:length(D));
    nrted=zerocross(Daux);
    nrted=mpicd+nrted-1;
    Daux=flipud(D(1:mpici));
    nrtei=zerocross(Daux);
    nrtei=mpici-nrtei+1;
    prr=1.4;
    if (abs(D(mpicd))>prr*abs(D(mpici))|(PKni-nrtei)>(nrted-PKni))
        % ---- PKni corresponds to Q wave, R wave will be on the right
        % ----
        Daux=D(mpicd:length(D));
        ncero=zerocross(Daux);
        if isempty(ncero) return; end
        ncero=mpicd+ncero-1;
        Ir=find(I>ncero); 
         if isempty(Ir) return; end
        mpda=I(Ir(1));
        if (ncero-PKni)/Fs>150e-3|(-1)*D(mpda)<ydd/10
            Sgran=1;
        else Rp=ncero;
        end
    else Sgran=1;
    end
    if Sgran==1
        % ---- PKni corresponds to S wave, R wave will be on the left ----
        Daux=flipud(D(1:mpici));
        ncero=zerocross(Daux);
        if isempty(ncero) return; end
        ncero=mpici-ncero+1;
        ilim=ncero-round(60e-3*Fs);
        if ilim<=0 ilim=1; end
        Daux=D(ilim:ncero);
        if (~isempty(Daux))
            [ymax2,imax2]=max(Daux);
            imax2=ilim+imax2-1;
            %if (PKni-ncero)/Fs>140e-3|(ymax2<(-1)*ydi/10)|((Xpb(ncero))<-abs(Xpb(PKni))*1/6)
            if (PKni-ncero)/Fs>140e-3|(ymax2<(-1)*ydi/100)|(Xpb(ncero)<-abs(Xpb(PKni))*1/3)  			%JGM
                Rp=PKni;
                type=4;  % ---- QS type ----
            else Rp=ncero;
            end
        end
    end
end


% ---- R wave location in the normal QRS type ----
if (type==1)
    Rp=PKni;
    R2p=[];
end

