%function x = slit_fn(n_fns, n_pix, noise)
%cannot make this a function when length, centre and meas_fn are all seperate arrays

% slit_fn(n_fns, n_pix, noise)
%	This function takes some number (n_fns) of experimentally determined slit
%	functions (meas_fns) and combines them to create a total slit function for
%	our spectrometer.

%	ARGUMENTS:
%			n_fns			# of functions measured
%			n_pix				total # of pixels (usually 2000)
%			noise				noise threshold for signal (signal below this
%								value is considered to be zero)
%	VARIABLES:
%			meas_fns			matrix containing the measured slit functions
%								rows 		= measured function # (ie 1:n_fns)
%								columns 	= function value (1:s_len(row #))
%
%	INPUT FILES:
%			whatever files contain the slit function data
%	***

% use a structure to represent the measured slit functions:
%	struct meas_fn
%		length 		length of the function (ie span in pixel space)
%		centre		pixel # of the centre of the function (measured
%							from the first pixel)
%		index			index to data matrix



%		s_fn			values of the function (i.e. there are s_len of these)

% try a cell instead; structure doesn't work.
% CELL FORMAT:
%		rows index individual measured slit functions (i.e. 1:1:n_fns)
%		columns are as follows:
%				COLUMN		CONTENTS
%					1			length of function (ie span in pixel space)
%					2			pixel # of the centre of the function (measured
%									from the first pixel)
%					3			vector containing the function values (ie an array
%									of length 1:1:[length of function]

% OKAY. HERE IS THE FINAL ABSOLUTELY LAST DEFINATELY THE LAST ATTEMPT TO GET
% A USEFUL DATA STRUCTURE FOR THIS THING. Use a struct to store length and
% center and matrix index to a seperate matrix that stores the actual function
% data. A thousand curses upon MatLab.

% INITIALIZE:
n_fns = 3;
coord(i) = [4 7 12];
n_pix = 2000;
noise = 1;

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

% READ IN SLIT FNS. temp fix:
k = 1:10;
temp(k) = 3;
meas_fn(i,k) = 3;
length(i) = 10;
centre(i) = 2;
%coord(i) = [35 67 109 499 700 1012 1024 1300 1500 1984]

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

% cell stuff:
%meas_fn(i,3) = {temp};
%meas_fn(i,1) = {10};

%meas_fn(i,2) = {4};
%meas_fn(i).s_fn(1:10) = temp;
%meas_fn(i).s_len = 10;

% DEFINE NOISE THRESHOLD (AND SUBTRACT NOISE?)
meas_fn(i,k) = meas_fn(i,k) - noise;
%meas_fn(i,3) = meas_fn(i,3) - noise;
%meas_fn(i).s_fn(1:meas_fn(i).s_len) = meas_fn(i).s_fn(1:meas_fn(i).s_len) - noise;
%s_fns(i,1:s_ln(i)) = meas_fns(i,1:s_ln(i)) - 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:length(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


% INTERPOLATE TO GET A FUNCTION FOR EACH PIXEL, AND OUTPUT THE RESULT
% (MATRIX OF PIXELS x FUNCTION VALUES)

% 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);

%temp:
n_pix = 20;

max_len = max(length(i));
s_fn = zeros(n_pix, max_len);

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"
