function [results] = get_fitting_concs_RC(SOI, results, core, D, z_pc, thresh, to_save)

% this function originally calculated RMSE as a measure of misfit. However in this revised version
% it uses a reduced chi calculation (98%). As both cores have 4 measurements there are 3 degrees of freedom It adds results to
% the struct intially produced in all_model_concs.m and saves as a new
% mat.file if to_save == 1.

% check core being analysed
if core == 2
    name = 'BH02';
elseif core == 3 
    name = 'BH03';
else 
    name = 'BH06';
end

d=D.(name); % set name
numits = length(results.N14_final); % get numits


% get predictions of concs at sample mid-depths for samples of interest.
mid_pred = (SOI * 3) - 1; % defines sequence of mid-depth predictions within results [2, 5, 8, 11 and so on if there were more samples]
mid_predictions = results.N14_final(:, [mid_pred]);

% how many degrees of freedom (revision from paper)
dof = length(SOI) - 1;

% The preprint used RMSE as a simple measure of misfit between model (mid_predictions and
% data *d.n14_concs. This isn't error weighted but with only a few
% measurements and reasonably consistent errors this is probably OK.

%         M2 = zeros(length(numits),1);
%         for i = 1:numits
%         M2(i) = sqrt(mean(mid_predictions(i,:) - d.n14_concs([SOI])') .^2);
%         end
%         M2 = M2';

    % These lines now do a reduced chi squared misfit
        M2 = zeros(length(numits),1);
            for j = 1:numits
            M2(j) = sum(((mid_predictions(j,:)-d.n14_concs([SOI])')./ d.dN14([SOI])').^2);
            end
        M2 = (M2./dof)';

        %  r chi from table C4 in Bevington and Robinson (a bit static if
        %  you ever want to assess a core with more than 3/4 samples!)
        if dof == 3
        Rchi95 = 2.605;
        Rchi98 = 3.279;
        elseif dof == 2
        Rchi95 = 2.996;
        Rchi98 = 3.912;
        end
        
        if thresh == 95
        Mrchi = M2(M2<=Rchi95);
        
        else
        Mrchi = M2(M2<=Rchi98);
        end

        
% get iterations in ascending order of fit to return
[~, minpc] = sort(M2,'ascend');
% np = n/100;
% ind = numits * np;
minpc2 = minpc(1:length(Mrchi)); % this is now indices (in orignal results) of best fitting n% of iterations

results.N14_clip10 = results.N14_final(minpc2,:); % N.B results are in order of fit which is handy
results.Ts_clip10 = results.Ts(minpc2,:); %in order of misfit
results.Zs_clip10 = results.Zs(minpc2,:); %in order of misfit
results.core = name;
results.numits = numits;
% results.bounds = bounds;
% results.M = M;
results.M2 = M2;
results.Mclip = M2(minpc2,:); % misfit in order
results.date = datestr(now, 'dd-mmm-yyyy'); 

% now save if that option selected
if to_save == 1
save_name = strcat('fitting_model_results','_',name,'_',datestr(now, 'dd-mmm-yyyy'),'_',num2str(numits),'_',num2str(thresh),'_fix_RChi');
save(save_name, 'results')
 
end

end