function focus_bootstrap4(filename)
% Focus by fitting to find the x corresponding to
% minimum spot size.
%
% The input data file must be in Octave text format (see data.txt).
%
% Based on focus_bootstrap3.m.
% Updated and cleaned up the code. Now running under Octave 3.
%     --- CCK 2009-Apr-29
%
% 2016-Apr-15 Switched to polyder instead of obsolete polyder. Swatted one (last?) bug.

degree = 2 % polynomial degree for the fit (typically 2).
N_trials = 10000;
confidence_level = 0.98;

% global mean_square_radius position   % share data with merit function 
%  (not needed for this version of the code)

% Read in data file
load('-text', filename, 'data')
position = data(:,1);
mean_square_radius = data(:,2);
N_points = prod(size(position))

figure(1)
hold off
plot(position, sqrt(mean_square_radius),'ro'); % plot the rms radii.
title('Focus Data')
xlabel('focuser position (in)')
ylabel('RMS spot radius (pixels)');
% print([filename, '_FocusData.pdf'])

figure(2)
hold off
plot(position,mean_square_radius,'ro'); %Plot the data
title('Best Fit')
xlabel('focuser position (in)')
ylabel('mean square spot radius')
hold on

% Find best fit polynomial from noisy data.
coeffs = polyfit(position, mean_square_radius, degree)';
   % Note that the transpose is used since coeffs returns a column vector.
% Take the derivative of the best fit polynomial.
deriv = polyder(coeffs);

% Generate best fit curve, plot for comparison
xmin = min(position);
xmax = max(position);
dx = (xmax - xmin)/100;
x = xmin:dx:xmax;
fit_I = polyval(coeffs, x);
plot(x,fit_I,'b-'); % Plot the fit
legend('data','best fit') % Add a legend, now that the plot is complete.
legend('location','southeast') % Better place for the legend
legend('boxon') % I like a box around the legend.
% print([filename, '_BestFit.pdf'])

% X0, where P(x) is minimized for xmin < x < xmax, is best focus position.
all_roots = roots(deriv) 
real_roots = all_roots(find(all_roots == real(all_roots)))
roots_lt_xmax = real_roots(find(real_roots<xmax));
roots_in_domain = roots_lt_xmax(find(roots_lt_xmax>xmin))
[Imax, i] = min(polyval(coeffs,roots_in_domain)); 
   % CCK 20050713 changed from max to min
X0 = roots_in_domain(i)



% ************** And now for the bootstrap! **********************
t_start = time;
start_bootstrap = ctime(time)

% Save the raw data. Must overwrite these arrays with randomly drawn
% data so badness.m can use them.
position_orig = position;
mean_square_radius_orig = mean_square_radius;

for iteration = 1:N_trials
   % pick N_points points, with replacement
   indeces = ceil( N_points*rand(N_points,1) );
   position  =  position_orig(indeces);
   mean_square_radius = mean_square_radius_orig(indeces);

   % do the fit

   % Find best fit polynomial from noisy data.
   coeffs = polyfit(position, mean_square_radius, degree)';
      % Note that the transpose is used since coeffs returns a column
      % vector (contrary to the octave documentation!).
   % Take the derivative of the best fit polynomial.
   deriv = polyder(coeffs);

   % X0, where P(x) is min for xmin < x < xmax, is best focus position.
   all_roots = roots(deriv) ;
   real_roots = all_roots(find(all_roots == real(all_roots)));
   roots_lt_xmax = real_roots(find(real_roots<xmax));
   roots_in_domain = roots_lt_xmax(find(roots_lt_xmax>xmin));
   [Imax, i] = min(polyval(coeffs,roots_in_domain));
   X0 = roots_in_domain(i);
   % CCK 2016-4-15 changed max to min in above block (think?)

   xs(iteration) = X0; % record the fit
   
end

t_finish = time;
bootstrap_completed = ctime(t_finish)
bootstrap_time_secs = t_finish - t_start
% bootstrap_time_mins = bootstrap_time_secs/60
% bootstrap_time_hrs  = bootstrap_time_mins/60


% Statistics
N_trials
mean_position = mean(xs)
standard_deviation = std(xs)
median_position = median(xs)

% Work out the confidence interval
confidence_level
xs = sort(xs);
lo_index = floor(0.5*(1-confidence_level)*N_trials);
hi_index = ceil(N_trials - 1 - 0.5*(1-confidence_level)*N_trials);
confidence_interval = [xs(lo_index), xs(hi_index)]
minimum_fit = min(xs)
maximum_fit = max(xs)

save -text bootstrap_log.txt
hold off
%end of main program




