% slit_fn
%	This script takes experimentally determined slit functions (meas_fns) 
%	and combines them to create a total slit function for our spectrometer.
%
%	For 2000 pixels and 4 useful slit functions, the script took 2379 seconds
%	(40 minutes) to execute. The final results, a matrix containing the slit
%	functions, is about 2 megabytes for functions that are around 120 pixels
%	wide. (July 7th)

%	LAST UPDATE: July 5th, 00

%	VARIABLES:
%	Input variables:
%		Read from slit function data files:
%	(g)	coord			array giving the first pixel # (ie lowest pixel #)
%								of each slit funtion
%	(g)	len			array containing the width (in pixels) of each slit function
%	(g)	meas_fn		matrix containing the measured slit functions
%								rows 		= measured function # (ie 1:n_fns)
%								columns 	= function value (1:length(row #))
%		Other input variables:
%			wv_fn			polynomial giving pixel # as function of wavelength
%			n_fns			# of measured slit functions
%	(g)	n_pix			# of pixels on the CCD
%			dc_noise		CCD dark current noise
%			bias_noise	CCD bias noise
%	(g)	h_width		desired half width of slit functions
%	(g)	max_iterations		max # of iterations to do before quitting
%			slit_fn_op				structure storing flags for various options:
%				slit_fn_op.norm		1 = normalize, 0 = don't
%				slit_fn_op.noise		1 = subtract dc and bias noise, 0 = don't
%				slit_fn_op.mf_disp	1 = display plot of measure slit functions, 0 = don't
%				slit_fn_op.compare	1 = read in another measured list function to compare
%								with interpolated ones, 0 = don't
%				slit_fn_op.stopwatch	1 = time the script, 0 = don't
%	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. 
%			max_cent 	centre coord of measured slit function closest to n_pix
%			min_cent		centre coord of measured slit function closest to 1
%	(g)	int_time		integration time of spectra
%	(g)	mean_TCCD	mean temperature of CCD while spectra was taken
%
%	--> 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 (.csv files)
%			csv_lis.lis			list of filenames of the data files (see the
%									m-file read_csv_slit_fns.m for more details)
%			bias_lookup.mat	generated by bias_fit.m, loaded by bias.m
%			dc_lookup.mat		generated by dc_fit.m, loaded by darkcurrent.m
%	CALLS:
%		find_centre				finds optimal centre pixel for each slit function
%		read_csv_slit_fns		reads in slit function data from .CSV files
%									(these are spectrometer raw output)
%		darkcurrent				calculates dark current noise for each pixel at a 
%									given mean CCD temperature
%		bias						calculates bias noise for each pixel at a 
%									given mean CCD temperature
%	SUMMARY OF FILES REQUIRED:
%		slit_fn.m		(driving program)
%		read_csv_slit.fns.m
%		find_centre.m
%		darkcurrent.m
%		bias.m
%		bias_fit.m
%		dc_fit.m
%		bias_lookup.mat
%		dc_lookup.mat
%
%	And for convolution of the spectrum, use spec_gen.m

% INITIALIZE:
clear n_fns n_pix dc_noise bias_noise i coord centre meas_fn s_fn
clear max_len lhs_len rhs_len min_temp s_loc len
clear z1 z2 a1 a2 m n 

clear global len meas_fn coord 

global n_pix meas_fn len h_width max_iterations coord int_time %mean_TCCD
global slit_fn_op		% program options, see below

% Create global struct to store slitfunction filenames, and the number of 
% seperate spectra in each of them that are averaged to get a slit function:
global sf_files


% Initialize options:
slit_fn_op.norm 		= 1;
slit_fn_op.noise 		= 1;
slit_fn_op.mf_disp 	= 1;
slit_fn_op.compare	= 0;
slit_fn_op.stopwatch = 1;

slit_fn_op;		% display options to user

if slit_fn_op.stopwatch 
   tic 
end

n_pix = 2000;

h_width = 50;	% desired half width, in pixels, of the slit functions
max_iterations = 1;

% Read in slit function data from CSV files
dir_path = 'c:\james\analyse\slit\slit_data\';
%dir_path = 'c:\james\analyse\slit\slit_data\do_not_use\';

wv_fn = [-8.5631 4884];		% polynomial for converting wave centre to pixel centre
									% note: this poly is NOT the same thing as the usual
                           % p = a + b*lambda relationship (in which a < 0 and
                           % b > 0).
n_fns = read_csv_slit_fns(dir_path, wv_fn, h_width);	% see READ_CSV_SLIT_FNS.M
i = 1:1:n_fns;
if n_fns < 1
   disp('SLIT_FN.M: Error, problem reading in slit functions.')
   return
else
   disp(strcat('SLIT_FN.M: No. of measured slit functions = ', num2str(n_fns)))
end

save noisy meas_fn
%return		% debugging

% Bring the noise: --> shift this to the read_csv_slit_fns2 procedure
%if slit_fn_op.noise
%   for m = 1:1:n_fns
%      dc_noise(m,:) = darkcurrent(mean_TCCD(m), coord(m), coord(m) + len(m) - 1);
%      bias_noise(m,:) = bias(mean_TCCD(m), coord(m), coord(m) + len(m) - 1);
%      k = 1:1:len(m);
%      meas_fn(m,k) = meas_fn(m,k) - dc_noise(m,k) - bias_noise(m,k);
%      meas_fn(m,k) = max(meas_fn(m,k), 0);		% cast out thy negative values
%   end	% for m = 1:1:n_fns
%end		% if slit_fn_op.noise

centre = coord;	% initialize centre to values that are definitely not the centre
lower(1:n_fns) = 1;

% Find centre of each slit function:
for m = 1:1:n_fns
	centre(m) = find_centre(m, centre(m), lower(m), len(m), 0);
end

%disp('printing from slit_fn.m, centre is...')
%centre

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


% get max slit function width
max_len = max(len);

% 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);					% remember, centre = RELATIVE centre (in pixels) 
												%           coord = ABSOLUTE centre (in pixels)
rhs_len = max(max_len - centre);
shift = lhs_len - centre;
centre = centre + shift;
clear temp
temp = meas_fn;
for m = 1:1:n_fns
   for k = 1:1:len(m)
      if (k + shift(m) > 0)	%bandaid fix
   		meas_fn(m,k + shift(m)) = temp(m,k);
      end
   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                  
                     
% Update centre and max_len:
centre = centre(1);
max_len = centre + rhs_len;
disp(strcat('Updated centre and max_len: centre = ', num2str(centre), ...
   ', max_len = ', num2str(max_len)));                     


% Now sort meas_fn so that the measured slit functions are organized in order
% of increasing absolute pixel centre:
disp(['Sorting slit functions in order of increasing pixel centre.'])
clear temp
p(i) = zeros(1,n_fns);
s_loc(i) = zeros(1, n_fns);
temp = coord;	
max_temp = max(temp);
% Check that no two values of temp(i) are the same:
	% MatLab note: you cannot change the upper limit on a for loop while the
	% loop is executing.
k = 1;
while k <= n_fns
%for k = 1:1:n_fns
   m = k + 1;
   while m <= n_fns
   %for m = (k+1):1:n_fns
      if temp(k) == temp(m)
         temp(m) = [];
         meas_fn(m,:) = [];		% delete redundant row of meas_fn 
         len(m) = [];				% and of len
         % Update n_fns and i:
         n_fns = n_fns - 1;
         %disp(strcat('k = ',num2str(k),', m = ',num2str(m)))  % debugging
         i = [1:n_fns];
      end
   	m = m + 1;
   end	% end of "while m <= n_fns"
	k = k + 1;   
end
if n_fns < length(coord)	% if n_fns has changed
   disp('Some of the measured slit functions were redundant and have')
   disp(strcat('been discarded. The new number of functions is n_fns = ',num2str(n_fns)))
   disp('	')
end
coord = temp;	% Update coord
% Sort slit functions in increasing order:
clear p
for k = 1:1:n_fns
   [p(k), s_loc(k)] = min(temp);
   temp(s_loc(k)) = max_temp + 1;   
end
clear temp
temp = meas_fn;
k = 1:1:n_fns;
meas_fn(k,1:max_len) = temp(s_loc(k),1:max_len);
%for k = 1:1:n_fns
%   meas_fn(k,1:max_len) = temp(s_loc(k),1:max_len);
%end
coord = p;
clear temp p

% Now we have coord and meas_fn as they were, except that the function indices
% (the rows of meas_fn) now increase along with the absolute centre pixel 
% coordinates of the functions. That is, meas_fn(1,1:max_len) gives the 
% measured slit function with lowest centre coordinate, meas_fn(2,1:max_len) 
% is the one with next lowest centre coordinate, and so on.

if slit_fn_op.norm
   method = 1;			% method = 1 	-> normalize by area
   						% 			= 2	->	normalize by height
   if method == 2
      % Normalize the slit functions by area:
      disp(['Normalizing slit functions by area.'])
		norm_area = (sum(meas_fn(1,:),2) + sum(meas_fn(2,:),2))/2;
		% arbitrarily choose 1st slit function as the standard
		j = 1:1:max_len;
		%k = 2:1:n_fns
		%clear temp_sum
		%temp_sum = sum(meas_fn(k,:),2);
		%meas_fn(k,j) = (norm_area/temp_sum).*meas_fn(k,j);
		% Vectorization of this loop is above
		for k = 2:1:n_fns
			clear temp_sum
		   temp_sum = sum(meas_fn(k,:),2);
		   meas_fn(k,j) = (norm_area/temp_sum).*meas_fn(k,j);
		end 
   	clear j temp_sum   
   else	
      % Normalize the slit functions by height:
      disp(['Normalizing slit functions by height.'])
      z(i) = max(meas_fn(i,:)');		% note transpose operator
      [norm_height  max_h] = max(z);
     	j = 1:1:max_len;
      for k = 1:1:n_fns
         scale_fac = norm_height/z(k);  
         meas_fn(k,j) = scale_fac.*meas_fn(k,j);
      end
   end
end	% if slit_fn_op.norm   

clear x

if slit_fn_op.mf_disp
   % display slit functions, probably for debugging purposes.
   disp(['Displaying measured slit functions.'])
	figure(3)
	clf
	hold on
	clear x
	for k = 1:1:n_fns
		x(k,:) = coord(k) - (centre - [1:max_len]);
	   plot(x(k,:), meas_fn(k,1:max_len),'r');
	end
	%x2 = coord(2) - (centre - [1:max_len]);
	%x3 = coord(3) - (centre - [1:max_len]);
	%plot(x1, meas_fn(1,1:max_len),'r');
	%plot(x2, meas_fn(2,1:max_len),'g');
	%plot(x3, meas_fn(3,1:max_len),'b');
	hold off
end	% if slit_fn_op.mf_disp

% DEBUGGING:
return

% Do interpolation to get the slit functions for each pixel:                       
method = 'linear';
disp(strcat('Degree of interpolation between measured slit functions = ', method));
n_pix = 280;
% initialize output slit function array (not necessary in MatLab, but speeds up the code)
s_fn = zeros(n_pix, max_len);



%n_pix = 241;	% speed up code; I only want a few interpolated functions right now (July 7)
n_pix = 2000;

%for j = 141:20:n_pix
for j = 1:1:n_pix
	for k = 1:1:max_len
      clear s_temp p_temp
      %s_temp(i) = meas_fn(i,k);		% vectorized code
      %p_temp(i) = coord(i) - (centre - k);
      for m = 1:1:n_fns
         s_temp(m) = meas_fn(m,k);
         p_temp(m) = coord(m) - (centre - k);
      end
      s_fn(j,k) = interp1(p_temp, s_temp, j - (centre - k), method);
   end
end

save slitfunctions s_fn max_len n_pix centre coord meas_fn n_fns %filename
% s_fn, max_len, n_pix, coord and centre are the variables required 
% to generate a convolved spectrum (see SPEC_GEN.M).
%
% meas_fn and n_fns just tell us about the measured slit 
% functions (it is probably useful to save them because
% at some later point, we might want to know what the measured
% functions were).

%Plot 2 measured fns:
%m = coord(1);	% coord of centre
%plot([(m - centre + 1):(m - centre + max_len)], s_fn(m,1:max_len), 'r');
%m = coord(2);
%plot([(m - centre + 1):(m - centre + max_len)], s_fn(m,1:max_len), 'b');

%m = abs(coord(2) - coord(1))/2;
%m = m - mod(m,1);
%m = min(coord(1),coord(2)) + m + 1;

if slit_fn_op.compare
   % plot some results:
	figure(1);
	clf;
	hold on;
   % compare a measured fn to an interpolated one:
	z1 = wk1read('s_230.wk1');
	z2 = z1';
	len(n_fns + 1) = length(z2);
	k = 1:1:len(n_fns + 1);
	noise = 1450;
	z2(k) = z2(k) - noise;
	meas_fn(n_fns + 1, k) = z2(k);
	cent_z2 = 230;
	cent_z2 = find_centre(n_fns + 1, cent_z2, 1, len(n_fns + 1), 0);
	coord(n_fns + 1) = 230 + cent_z2;
	m = coord(n_fns + 1);
	plot([(m - cent_z2 + 1):(m - cent_z2 + len(n_fns+1))], z2(1:len(n_fns+1)), 'g');

	% Plot interpolated function(s):
	for m = 141:20:241
		%m = 141;
		n = m;
		plot([(n - centre + 1):(n - centre + max_len)], s_fn(m,1:max_len), 'k');
	end

	title('BLACK = interpolated function, GREEN = measured function');
end	% if slit_fn_op.compare

if slit_fn_op.stopwatch
   format long
	toc
	format short
end

%n = coord(2);
%plot([(n - centre + 1):(n - centre + max_len)], s_fn(m,1:max_len), 'g');
%z1 = wk1read(strcat(dir_path,'s_230.wk1'));




