% slit_fn
%	This script takes experimentally determined slit functions (meas_fns) 
%	and combines them to create a total slit function for our spectrometer.

%	VARIABLES:
%	Input variables:
%		Read from slit function data files:
%			coord			array giving the first pixel # (ie lowest pixel #)
%								of each slit funtion
%			length		array containing the width (in pixels) of each slit function
%			meas_fn		matrix containing the measured slit functions
%								rows 		= measured function # (ie 1:n_fns)
%								columns 	= function value (1:length(row #))
%		Other input variables:
%			n_fns			# of measured slit functions
%			n_pix			# of pixels on the CCD
%			noise			# noise threshold
%
%	Output and other variables
%			centre		array giving the centre of "mass" of each slit function. 
%								After the script has run, this will be the same for 
%								all pixels.
%			coord			Okay, this is the same variable as mentioned above. But
%								after the script has run, it contains the pixel 
%								coordinate of the centre of each measured slit function.
%			s_fn			output matrix for the slit functions
%								rows 		= pixel #
%								columns 	= function value (1:max_len)
%			max_len		maximum width of measured slit functions. 
%
%	--> So, just to summarize, the key output of this script is contained in the 
%	variables "s_fn", "centre" and "max_len".
%
%	INPUT FILES:
%			whatever files contain the slit function data

% INITIALIZE:
n_fns = 3;
n_pix = 2000;
noise = 1;

%temp:
n_pix = 2000;

i = 1:1:n_fns;
%j = 1:1:n_pix;

coord(i) = [4 7 12];

% READ IN SLIT FNS. temp fix:
k = 1:10;
temp(k) = 3;
meas_fn(i,k) = 3;
len(i) = 10;
centre(i) = 2;

% note: use the fn "wk1read('filename')" to read in Lotus (wk1) spreadsheet data

meas_fn(1,k) = 3;
meas_fn(2,k) = 3;
meas_fn(1,1) = 10;
meas_fn(2,8) = 10;

% DEFINE NOISE THRESHOLD (AND SUBTRACT NOISE?)
meas_fn(i,k) = meas_fn(i,k) - noise;

% FIND WEIGHTED CENTRES (ie centre of "mass" of the function)
clear temp temp_k k centre
temp_k(i) = 0;
temp(i) = 0;
for m = 1:1:n_fns
   for k = 1:1:len(m)
   	temp_k(m) = temp_k(m) + k.*meas_fn(m,k);
   	temp(m) = temp(m) + meas_fn(m,k);
   end
end
centre(i) = temp_k(i)./temp(i);

clear temp
temp(i) = mod(centre(i),1);
% Round off r_cm to get pixel numbers:
for m = 1:1:n_fns
	if temp(m) < 0.5
		centre(m) = centre(m) - temp(m);		 	% round down
	else
   	centre(m) = centre(m) + 1.0 - temp(m);	% round up
	end
end

% Shift coords of measured functions so that coord now stores the pixel
% position of the centre of each one:
coord(i) = coord(i) + centre(i);
% note, the above line MUST be executed before the shifting code below runs!!

% get max slit function width
max_len = max(len(i));
% initialize output slit function array (not necessary in MatLab, but speeds up the code)
s_fn = zeros(n_pix, max_len);

meas_fn	%debugging output

% To account for the possibility of measured slit functions being of different
% pixel widths, we should find the maximum function width and then shift the 
% other functions so that they are centred within that pixel range. That is,
% we shift the functions so that they are all centred at the same column index
% of meas_fn:
lhs_len = max(centre);
rhs_len = max(max_len - centre);
shift(i) = lhs_len - centre(i);
centre(i) = centre(i) + shift(i);
clear temp
temp = meas_fn;
for m = 1:1:n_fns
	for k = 1:1:len(m)
   	meas_fn(m,k + shift(m)) = temp(m,k);
	end
   for k = 1:1:shift(m)
      meas_fn(i,k) = 0;
   end
   for k = (len(i) + shift(i) + 1):1:(centre(i) + rhs_len)
      meas_fn(i,k) = 0;
   end
end						% note to self: this code redefines array lengths a lot,
                     % and so it's probably not as fast as it could be if the
                     % array lengths were preallocated in memory. FIX LATER
                     
                     
% Do linear interpolation to get the slit functions for each pixel:                     
for j = 1:1:n_pix
   % find the two closest measured functions to pixel j:
   clear temp
   temp(i) = abs(coord(i) - j);
   [min_temp(1), s_loc(1)] = min(temp);
   temp(s_loc(1)) = max(temp);
   [min_temp(2), s_loc(2)] = min(temp);

   if ((j < coord(s_loc(1))) & (j < coord(s_loc(2)))) | ((j > coord(s_loc(1))) & (j > coord(s_loc(2)))) | (min_temp(1) == 0)
      s_fn(j,1:max_len) = meas_fn(s_loc(1),1:max_len);
   else
      % interpolate to get the function values
      if (j < coord(s_loc(1)))
         x1 = coord(s_loc(1)) - j;
         x2 = j - coord(s_loc(2));
      else
         x2 = coord(s_loc(2)) - j;
         x1 = j - coord(s_loc(1));
      end
      for m = 1:1:max_len
      	s_fn(j,m) = (meas_fn(s_loc(1),m)*x2 + meas_fn(s_loc(2),m)*x1)/(x1 + x2);
      end
   end 	
end % end of "for j = 1:1:n_pix"

centre = centre(1);
