function x = find_centre(m, centre, lower, upper, it_count);
% find optimal weighted centres for one slit function, indexed by m. Recursive.

%	ARGUMENTS:
%	note: all arguments are NUMBERS, not vectors or matrices
%		m				index of the slit function being centred (between 1 and n_fns)
%		centre		current centre of slit function
%		lower			current lowermost pixel of slit function
%		upper			current uppermost pixel of slit function
%				NOTE: "lower" and "upper" are both relative pixel coords,
%				absolute pixel coords (ie they are not pixel # on the CCD,
%				but indices of pixels within "meas_fn")
%		it_count		keeps track of recursive calls to this function. We want
%						to stop iterating after reaching max_iterations calls (just
%						because it is assumed that we're not making progress after
%						that point).
%	Global variables required:
%		meas_fn		measured slit function data, stored in an n_fns x len(m) matrix
%		len			pixel width of measured slit functions
%		h_width		desired half width of slit functions
%		max_iterations		max # of iterations to do before quitting
%
%	CALLS:			none
%	CALLED BY:		slit_fn.m

global meas_fn
global len
global h_width
global max_iterations

clear temp temp_k k

if it_count == max_iterations		% give up trying to optimize centre
   disp(strcat('Max # of iterations reached in FIND_CENTRE, it_count = ', num2str(max_iterations)))
   x = centre;
   return
end

temp_k = 0;
temp = 0;

for k = lower:1:upper
  	temp_k = temp_k + k.*meas_fn(m,k);
  	temp = temp + meas_fn(m,k);
end
z = temp_k/temp;	% z = our new centre coordinate
% Round off the new centre:
if (mod(z,1) < 0.5)
   z = z - mod(z,1);
else
   z = z + 1 - mod(z,1);
end

if (z == centre)	% then we are finished
   x = z;
else	% ie if (z ~= centre), we are not finished
   lower = z - h_width;
   upper = z + h_width;
   % Round off to integer pixel numbers:
   if (mod(lower,1) < 0.5)
      lower = lower - mod(lower,1);
   else
      lower = lower + 1 - mod(lower,1);
   end
   if (mod(upper,1) < 0.5)
      upper = upper - mod(upper,1);
   else
      upper = upper + 1 - mod(upper,1);
   end
   % Make sure that lower(m) and upper(m) are in range:
   lower = max(1, lower);
   upper = min(len(m), upper);
   % Look for new centre:
   x = find_centre(m, z, lower, upper, it_count + 1);    
end
