function D3M_3Dfun(RR,videoname, foldername, T, Tshift, step)
% D3M-3Dfun(RR) - creates 2D contour & 3D surface D-3M movie frames
%
% input --> RR: two-column matrix containing
%           in the first column the time instants of the R peaks (in s)
%           in the second column the RR intervals (in s)
%           videoname: string containing the name for the 
%           video
%           foldername (optional): string containing the name of the
%           folder where the frames will be saved
%           T (default=10): time interval to display in one frame (in min)
%           Tshift (default=1): time shift between frames (in min)
%           step (default=0.00015): step between plot contour lines
%
% Developed by A Burykin (burykin@gmail.com), L Citi, 
% MD Costa and AL Goldberger
% and modified by S Mariani (sara.mariani@wyss.harvard.edu)
% Wyss Institute at Harvard
% when using this function, please reference: Burykin, Anton et al.
% "Dynamical density delay maps: simple, new method for visualising 
% the behaviour of complex systems." BMC medical informatics and 
% decision making 14.1 (2014): 6.

narginchk(2, 6);

if (nargin>=3 && ~isempty (foldername))
    % save frames in folder specified by user
    savepic=1;
    mkdir(foldername);       
    else savepic=0;
end

if (nargin<4 || isempty (T))
    % use default interval
    T=10;
end
    
if (nargin<5 || isempty (Tshift))
    % use default shift
    Tshift=1;
end

if (nargin<6 || isempty (step))
    % use default step
    step=0.00015;
end

vidObj = VideoWriter(videoname);
vidObj.FrameRate = 10;
open(vidObj);


t=RR(:,1);
t=t-t(1);
t=t/60.0; % sec2min
rr=RR(:,2);

scrsz = get(0,'ScreenSize');

k1=1;
k2=find(t>=T,1);

lolim=min(RR(:,2));
hilim=max(RR(:,2));

k=1;

while t(k2)<=t(end)

    t1=t(k1);
    t2=t(k2);

    t_disp=t(k1:k2);
    rr_disp=rr(k1:k2);

    y1=rr_disp(1:end-1);
    y2=rr_disp(2:end);

    figure
    [ctrs1, ctrs2, F, hAxes] = dscatter2(y1,y2,'PLOTTYPE','contour');
    close

    fig1=figure('Position',...
    [0.05*scrsz(3) 0.05*scrsz(4) 0.8*scrsz(3) 0.89*scrsz(4)],...
    'Color',[0 0 0]);
    colormap jet
       
    subplot(4,1,[1 2])
    surf(ctrs1, ctrs2, F, 'FaceColor','interp',...
    'EdgeColor','none',...
    'FaceLighting','phong')
    xlim([0.5 1.4])
    ylim([0.5 1.4])
    colormap jet
    view(5,70)
    axis off
    subplot(4,1, 3)
    contour(ctrs1, ctrs2, F,'LineWidth',2,'LevelStep',step)
    view(5,89)
    colormap jet

    xlim([0.5 1.4])
    ylim([0.5 1.4])
    zlim([0.0  0.001])
    axis off
   
    subplot(4,1,4)
    plot(t_disp,rr_disp,'r','LineWidth',2.5)
    xlim([(min(t_disp)-1) (max(t_disp)+2)])
    ylim([0.5 1.4])
    axis off
 
    % capture frame here  
    Fr1=getframe(fig1);
    Im1=frame2im(Fr1);
    
    if savepic
    im_file=[foldername '/' int2str(k) '.bmp'];
    imwrite(Im1,im_file);
    end

    writeVideo(vidObj,Fr1);
    clear Fr1 Im1 Im1_map im_file;
    close(fig1);
    
    k1=find(t>=t1+Tshift,1);
    k2=find(t>=t2+Tshift,1);
end
close(vidObj);