function  [spurt_data_subject graph_handle] = find_spurts(simulation_parameters, id)
    %
    % This function finds the spurts and optionally plots data for each curve
    %
    % Input parameters
    % ----------------
    % simulation_parameters: see description of simulation_parameters in
    % file global_constants.m
    % 
    % Output values
    % -------------
    % spurt_data: see description of spurt_data_template in file 
    % global_constants.m
    %
    % Note: stderr_location, stderr_amplitude and stderr_duration are set to -1
    % if compute_std_errors is set to false.
    %

    % load global constants
    global_constants;

    if (id > size(simulation_parameters.subject_ids,1))
        error(['No subject id data (simulation_parameters.subject_ids) found for subject id ', num2str(id)])
    end
    subject_id = simulation_parameters.subject_ids(id);

    disp(['Processing case id: ', num2str(subject_id)]);

    if (id > size(simulation_parameters.group_ids,1))
        error(['No group data (simulation_parameters.group_ids) found for subject id ', num2str(id)])
    end
    group = simulation_parameters.group_ids(id);  % retrieving group for this case
    
%     if (id > size(simulation_parameters.sample_ys,1))
%         error(['No data (simulation_parameters.sample_ys) data found for subject ', num2str(id)])
%     end
%     raw_data_sample = simulation_parameters.sample_ys(id,:);  % extract data for this case
    if (id > size(simulation_parameters.sample_ys))
        error(['No data (simulation_parameters.sample_ys) data found for subject ', num2str(id)])
    end
    raw_data_sample = simulation_parameters.sample_ys(:,id);  % extract data for this case
        
    if (size(simulation_parameters.sample_ts) == 0)
        error('No sample points (simulation_parameters.sample_ts) data found')
    end
    
    if (size(simulation_parameters.degree) == 0)
        error('degree not set')
    end

    if (size(simulation_parameters.pvalue) == 0)
        error('pvalue not set')
    end
    
    if (size(simulation_parameters.error_level) == 0)
        error('error_level not set')
    end

    if (size(simulation_parameters.output_all_maxima) == 0)
        error('output_all_maxima not set')
    end
    
    if (size(simulation_parameters.compute_std_errors) == 0)
        error('compute_std_errors not set')
    elseif (simulation_parameters.compute_std_errors == true)
        if (size(simulation_parameters.boot_samples) == 0)
            error('boot_samples value not set')
        end
    end
    
    if(size(raw_data_sample,2) ~= size(simulation_parameters.sample_ts,2))
        error('inconsistent number of samples in data and in sample_ts')
    end

    
    if (size(simulation_parameters.lambda_fct) == 0)
        error('lambda_fct not set')
    end
    
    if (size(simulation_parameters.lambda_stderr) == 0)
        error('lambda_stderr not set')
    end
    
    spurt_data = spurt_data_template;

    % fit smooth FDA curve to the raw_data. Also estimate standard error of
    % the fit at each raw x value
    [ hgtfd0 stderr y2cMap0 ] = create_fd_objects( simulation_parameters.sample_ts, ...
                                                   raw_data_sample, ...
                                                   simulation_parameters.lambda_fct, ...
                                                   simulation_parameters.lambda_stderr, ...
                                                   simulation_parameters.degree );

    % find the locations (i.e., x values) where maxima and minima in the
    % degree'th and degree+1'th derivatives of the fitted function (e.g., if
    % degree = 0, find max and mins in the function and the velocity)
    disp('Finding minima and maxima...');
    max_min_mat_loc = find_min_max_locations(simulation_parameters.sample_ts, hgtfd0, ...
        simulation_parameters.degree, simulation_parameters.error_level);
    x_values        = [max_min_mat_loc.location]';   % extract locations of the maxs and mins

    % compute the standard error for the degree'th and degree+1'th derivatives
    stderr_degree0 = get_stderr_estimates(x_values, hgtfd0, stderr, y2cMap0, ...
        simulation_parameters.degree);

    % determine which spurts are significant
    significant_spurt_entry = determine_spurts_significance(max_min_mat_loc, simulation_parameters.pvalue, ...
        hgtfd0, stderr_degree0, ...
        simulation_parameters.degree);

    % calculate spurt amplitude and duration
    disp('Estimating location, duration and amplitude...');

    if (simulation_parameters.compute_std_errors)
        % estimate std errors using the bootstrap method 
        [spurt_data_subject] = ...
            get_spurt_data_bootstrap( ...
                simulation_parameters.sample_ts, ...
                raw_data_sample,  ...
                significant_spurt_entry, ...
                hgtfd0, simulation_parameters.degree, ...
                simulation_parameters.boot_samples, ...
                simulation_parameters.lambda_fct, ...
                simulation_parameters.lambda_stderr, ...
                simulation_parameters.error_level, ...
                simulation_parameters.pvalue);

    else
        spurt_data_subject  = ...
            get_spurt_data_means_only(significant_spurt_entry, hgtfd0, ...
                simulation_parameters.degree, ...
                simulation_parameters.output_all_maxima);
    end
    
    for i=1:size(spurt_data_subject,2)
        spurt_data_subject(i).groupId          = group;
        spurt_data_subject(i).caseId           = simulation_parameters.subject_ids(id);
        spurt_data_subject(i).spurtId          = i;
    end
    
    if(size(simulation_parameters.group_labels) == 0)
        graph_title = ['Sample no. ', num2str(simulation_parameters.subject_ids(id))];
    else
        graph_title = [char(simulation_parameters.group_labels(group)), ' - Sample no. ', num2str(simulation_parameters.subject_ids(id))];
    end
    rng = [x_values(1), x_values(size(x_values,1))];  % range of x values covered in the raw data

    x_values = [spurt_data_subject.center_location];

    graph_handle = curve_display(hgtfd0, simulation_parameters.pvalue, rng, 0, stderr, y2cMap0, ...
        x_values, raw_data_sample, simulation_parameters.sample_ts, graph_title, ...
        simulation_parameters.x_label, ...
        simulation_parameters.y_label );

