function h = curve_display(data_fd, pvalue, rng, degree, stderr, y2cMap, ...
    max_t, y, t, graph_title, label_x, label_y)
    %
    % This function shows:
    %         function with the raw data
    %         residuals from of smoothed vs raw data 
    %         velocity     curves with 95% confidence bands
    %         acceleration curves with 95% confidence bands
    %         jerk curves with 95% confidence bands
    %         location of min and max in velocity and accelearation curves

    %  Mandatory input variables:
    %  data_fd  Functional data object of the smooth
    %  pvalue   Statistical significance to use
    %  rng      Range over which data_fd is defined 
    %  degree   Indicates which derivative of the data_fd corresponds to the
    %           base function. Generally, the base function is the 0th
    %           derivative, but in some cases it might be the first one, such
    %           as for height (growth), we might be interested in maxima of the
    %           velocity (rate of change), not maxima of the growth curve per se
    %  stderr   An FDA standard error matrix
    %  y2cMap   A matrix produced by SMOOTH_BASIS that maps
    %           a data vector into the corresponding coefficient
    %           vector
    %
    % Optional input variables:
    %  max_t  Vector of t values where significant maxima are located
    %  y   Vector of raw data (y values) to be plotted alongside with the function
    %  t   t values corresponding to data_mat
    %  graph_title    title of the graph
    %  label_x        string to use to label the x axis
    %  label_y        string to use to label the y axis

    % loading global constants
    global_constants;

    disp('Plotting for spurts and plateaus (including original function)...')

    h = figure('visible','off');

    max_count = size(max_t);
    if (max_count > 0)
        % now, get corresponding y values 
        fct_y  = eval_fd(max_t, data_fd, 0);
        vel_y  = eval_fd(max_t, data_fd, 1);
        acc_y  = eval_fd(max_t, data_fd, 2);
        jerk_y = eval_fd(max_t, data_fd, 3);
    else
        fct_y  = [];
        vel_y  = [];
        acc_y  = [];
        jerk_y = [];
    end

    %  get range of xvaluess and set up a fine set of values for plotting
    npoints   = 101;  % number of points for a smooth plot (dense enough)
    xvalueslo = rng(1);
    xvalueshi = rng(2);
    xvaluesfn = linspace(xvalueslo, xvalueshi, npoints)';  %  xvalues for plotting

    fd_basis  = getbasis(data_fd);

    %  evaluate data values, velocity, and acceleration for each 
    %  plotting value of xvalues
    %
    fcthat  = eval_fd(xvaluesfn, data_fd, degree);
    velhat  = eval_fd(xvaluesfn, data_fd, degree+1);
    acchat  = eval_fd(xvaluesfn, data_fd, degree+2);
    jerkhat = eval_fd(xvaluesfn, data_fd, degree+3);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    set(gcf, 'Color',[1,1,1])

    %  plot function values and raw data
    subplot(4,1,1)

    % computing confidence band
    parambasis  = eval_basis(xvaluesfn, fd_basis, degree);
    paramLmat   = parambasis*y2cMap;
    stderrparam = sqrt(diag(paramLmat * diag(stderr) * paramLmat'));
    fctband     = getBandWidthFactor(pvalue).*stderrparam;

    % adjust y display so that the displayed area is a bit larger than the
    % graph
    yvaluesLowConf  = fcthat - fctband;
    yvaluesHighConf = fcthat + fctband;
    yvalueslo       = min(yvaluesLowConf);
    yvalueshi       = max(yvaluesHighConf);
    
    if (size(y) > 0)
        plot(t, y, 'k.', ...
             xvaluesfn, yvaluesLowConf, 'k:', ...
             xvaluesfn, yvaluesHighConf, 'k:', ...
             xvaluesfn, fcthat, 'k-', ...
             max_t, fct_y, 'k*', ...
             max_t, fct_y, 'ko', ...
             rng, [0,0], 'k--')
    else
        plot(xvaluesfn, yvaluesLowConf, 'k:', ...
             xvaluesfn, yvaluesHighConf, 'k:', ...
             xvaluesfn, fcthat, 'k-', ...
             max_t, fct_y, 'k*', ...
             max_t, fct_y, 'ko', ...
             rng, [0,0], 'k--')
             
    end

    axis([xvalueslo,xvalueshi,yvalueslo,yvalueshi]);
    xlabel(label_x);
    ylabel(label_y);
    title(graph_title)


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % plotting velocity
    %
    subplot(4,1,2)

    % computing confidence band
    parambasis  = eval_basis(xvaluesfn, fd_basis, degree+1);
    paramLmat   = parambasis*y2cMap;
    stderrparam = sqrt(diag(paramLmat * diag(stderr) * paramLmat'));
    velband     = getBandWidthFactor(pvalue).*stderrparam;

    vellowbound  = velhat - velband;
    velhighbound = velhat + velband;

    if (size(y) > 0)
        % because data points are passed to the function, compute
        % and display corresponding points on the velocity curve
        deriv_vals = eval_fd(t, data_fd, degree+1);
        plot(t, deriv_vals, 'k.', ...
             xvaluesfn, velhat, 'k-', ...
             xvaluesfn, velhighbound, 'k:', ...
             xvaluesfn, vellowbound, 'k:', ...
             max_t, vel_y, 'k*', ...
             max_t, vel_y, 'ko', ...
             rng, [0,0], 'k--')
    else
        % Don't plot sampling value points
        plot(xvaluesfn, velhat, 'k-', ...
             xvaluesfn, velhighbound, 'k:', ...
             xvaluesfn, vellowbound, 'k:', ...
             max_t, vel_y, 'k*', ...
             max_t, vel_y, 'ko', ...
             rng, [0,0], 'k--')
    end

    yvello = min(vellowbound);
    yvelhi = max(velhighbound);
    axis([xvalueslo,xvalueshi,yvello,yvelhi])
    xlabel(label_x);
    ylabel('Velocity')

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % plot acceleration curve with confidence limits

    % computing confidence band
    parambasis  = eval_basis(xvaluesfn, fd_basis, degree+2);
    paramLmat   = parambasis*y2cMap;
    stderrparam = sqrt(diag(paramLmat * diag(stderr) * paramLmat'));
    accband     = getBandWidthFactor(pvalue).*stderrparam;

    subplot(4,1,3)
    acclowbound  = acchat - accband;
    acchighbound = acchat + accband;
    if (size(y) > 0)
        % because data points are passed to the function, compute
        % and display corresponding points on the acceleration curve
        acc_vals = eval_fd(t, data_fd, degree+2);
        plot(t, acc_vals, 'k.', ...
             xvaluesfn, acchat, 'k-', ...
             xvaluesfn, acchighbound, 'k:', ...
             xvaluesfn, acclowbound, 'k:', ...
             max_t, acc_y, 'k*', ...
             max_t, acc_y, 'ko', ...
             rng, [0,0], 'k--')
    else
        plot(xvaluesfn, acchat, 'k-', ...
             xvaluesfn, acchighbound, 'k:', ...
             xvaluesfn, acclowbound, 'k:', ...
             max_t, acc_y, 'k*', ...
             max_t, acc_y, 'ko', ...
             rng, [0,0], 'k--')
    end

    yacclo = min(acclowbound);
    yacchi = max(acchighbound);
    axis([xvalueslo,xvalueshi,yacclo,yacchi]),
    xlabel(label_x);
    ylabel('Acceleration')


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % plot jerk curve with confidence limits

    % computing confidence band
    parambasis  = eval_basis(xvaluesfn, fd_basis, degree+3);
    paramLmat   = parambasis*y2cMap;
    stderrparam = sqrt(diag(paramLmat * diag(stderr) * paramLmat'));
    jerkband    = getBandWidthFactor(pvalue).*stderrparam;

    subplot(4,1,4)
    jerklowbound  = jerkhat - jerkband;
    jerkhighbound = jerkhat + jerkband;

    if (size(y) > 0)
        % because data points are passed to the function, compute
        % and display corresponding points on the jerk curve
        jerk_vals = eval_fd(t, data_fd, degree+3);
        plot(t, jerk_vals, 'k.', ...
             xvaluesfn, jerkhat, 'k-', ...
             xvaluesfn, jerkhighbound, 'k:', ...
             xvaluesfn, jerklowbound, 'k:', ...
             max_t, jerk_y, 'k*', ...
             max_t, jerk_y, 'ko', ...
             rng, [0,0], 'k--')
    else
        plot(xvaluesfn, jerkhat, 'k-', ...
             xvaluesfn, jerkhighbound, 'k:', ...
             xvaluesfn, jerklowbound, 'k:', ...
             max_t, jerk_y, 'k*', ...
             max_t, jerk_y, 'ko', ...
             rng, [0,0], 'k--')
    end

    yjerklo = min(jerklowbound);
    yjerkhi = max(jerkhighbound);
    axis([xvalueslo,xvalueshi,yjerklo,yjerkhi]),
    xlabel(label_x);
    ylabel('Jerk')

    % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % % display phase plane plot of velocity and acceleration
    % subplot(4,2,2)
    % plot(velhat, acchat, '.k-')
    % xlabel('Velocity');  
    % ylabel('Acceleration');  
    % title('Phase Plane Plot')
    % 
    % if (size(t) > 0)
    %     subplot(4,2,4)
    %     fit_vals  = eval_fd(t, data_fd, degree);
    %     residuals = y - fit_vals;
    % 
    %     plot(t, residuals, 'k');
    %     xlabel(label_x);  
    %     ylabel('Residual');  
    %     title('Residuals (raw data - fit data)')
    %     
    %     subplot(4,2,6)
    %     hist(residuals);
    %     h = findobj(gca,'Type','patch');
    %     set(h,'FaceColor','k','EdgeColor','w')
    %     title('Residuals distribution')
    % end
    % 
%    disp('Press any key to continue')
%    pause
return;