function x = read_csv_slit_fns(dir_path, wv_fn)
% read_csv_slit_fns(dir_path)
%	This function reads in data from a series of .csv files (the 
%	raw data files produced by the spectrometer) to use for 
%	generating the total slit function for our spectrometer.

%	ARGUMENTS:
%		dir_path			path to directory containing the data files
%		wv_fn				polynomial giving the relation between the CCD's
%							wavelength centre and pixel centre (where does
%							it come from? There are some who call him...Kruger).
%
%	INPUT FILES:
%		csv_lis.lis		file containing a list with the filenames of csv
%							files to read.
%		ddmmyyhh.csv	data files (name gives date and time)
%	Required directory structure for input files:
%		DIRECTORY:				CONTENTS:
%		  slit					scripts: slit_fn.m, find_centre.m, read_csv_slit_fns.m
%			|						the list file, csv_lis.lis
%			|
%			|__[data 
%				directories]
%			
%	OUTPUT:
%		x = # of slit functions. Also sets the global variable coords to
%			contain the starting pixel of each slit function.
%		If the function fails, it returns x = -1. This can happen for the
%		following reasons:
%			- 	list of csv files is exhausted without finding any valid slit
%				functions
%
%	VARIABLES:
%			wv_cent		wavelength centre of slit fn
%			p_cent		pixel centre of slit fn (derived from wv_cent using wv_fn)
%	(g)	int_time		integration time of spectra
%			mean_TCCD	mean temperature of CCD while spectra was taken
%	(g)	n_pix			# of pixels on the CCD
%			T_min			min temperature for calculating dc and bias noise
%			T_max			max temperature for calculating dc and bias noise
%			wid			prescribed half width for the slit functions read in
%
%	CALLS:				darkcurrent.m
%							bias.m
%
%	CALLED BY:			slit_fn.m
%
%	COMMENTS:
%		- the script assumes 2000 pixels, ie there is no global n_pix variable
%		July 26 - modify so that functions are averaged over all spectra in the
%					 csv files, rather than using just the first one in the file.
%		- the value of wid MUST be less than 100!!!
%		--> July 27: stuff to add to this:
%DONE Jul28	- check CCD temp for each spectra, store them in global var
%				- maybe validate temps - check same within CSV file
%DONE Jul28	- subtract noise for each individual spectrum BEFORE
%					averaging over contents of each csv file
%DONE Jul30	- check for cutoff (saying, 15000 counts) and saturation (65540 counts)
%					in each individual spectrum, before averaging
%				- averaging - weight by counts?
%				- for peak normalization: by area should be good when all peaks
%					are averaged? Although, averaging over CSV files is kind of
%					arbitrary way to average. 
%						Instead, could perhaps store all spectra, from all CSV files,
%					individually and then average ones with ~same centre. 
%DONE Jul30	- does this code handle slit functions that fall off the beginning
%					or end of the ccd pixel range?
%DONE Jul30	- need to check centres of new functions against centres of functions
%					already read in, and add them to average of function already read
%					if they are the same.
%		--> what program does is to discard spectra in a given csv file if they 
%			don't have the same centre as spectra in that previous file. Jul30
%
% 	BUG		- 	for file 28060020.csv, something weird happens. Tried to fix, no
%					success, so just leave this file out for now. The functions
%					in it don't have the same centres anyways.

clear filename wv_cent
clear i n
clear meas_fn len coord

global n_pix meas_fn len coord int_time %mean_TCCD
%global dc_noise

global slit_fn_op		% program options
global sf_files		% struct for keeping track of slit fn filenames & # of 
							% averaged spectra that are used for each one
                     
sf_files.filename = [];		% records filenames of csv files that are used
									% (i.e. that provide acceptable slit lineshapes)                     
sf_files.n_avg = [];			% records # averages for given csv file
sf_files.temp = [];			% records mean CCD temperatures of spectra
sf_files.files_not_used = [];		% Useful for telling us which csv files to discard
                     
wid = h_width;		% Half width of slit functions to read in (just my guess
						% at a good number)
           			% *!*!* IMPORTANT *!*!*
            		% wid MUST be less than 100 for this program to run properly.
            
% Cutoff and saturation count levels:
counts_cutoff 	= 15000;
counts_sat 		= 65540;

% Max and min temperatures for which we can calculate the noise:            
T_max = 249;
T_min = 238;
disp('Temperature range for calculating noise:')
disp(['    ' 'Max:' '    ' 'Min:'])
disp(['    ' num2str(T_max) '      ' num2str(T_min)])


% open list file:
fid = fopen('csv_lis.lis','r');
if fid == -1
   disp('READ_CSV_SLIT_FN.M: Error, file CSV_LIS.LIS could not be opened.');
   x = -1; 
   return
end

i = 0;
while ~feof(fid)	% while not end of file
   i = i + 1;
   filename(i,1:8) = fscanf(fid,'%8s',1);
   filename(i,9:12) = '.csv';
   fgetl(fid);
end
fclose(fid);

n_fns = i;	% assume only one slit fn in each csv file. However, all of 
				% the spectra in each file are averaged to produce this one slit
            % function. It is assumed that all of the spectra in each csv
            % file have the same pixel centre. They are supposed to, so
            % this program doesn't check that they do.
				% 		Note, n_fns is a LOCAL variable - not the same as 
            % the one in SLIT_FN.M.
            
            
% Preallocate memory to speed up the code:
dc_noise 	= zeros(1,2*wid+1);
bias_noise 	= zeros(1,2*wid+1);
int_time 	= zeros(1,n_fns);
wv_cent 		= zeros(1,n_fns);
cent 			= zeros(1,n_fns);

%for i = 1:1:n_fns
i = 1;
p = 1; 	% counter for sf_files.files_not_used
disp(['Total no. of functions at start = ' num2str(n_fns)])
while i <= n_fns
   % Tell user what function we are reading:
   disp(['  '])
   disp(['Now reading function no. = ' num2str(i) ])
   
   % Open data file:
   fid = fopen(strcat(dir_path,filename(i,1:12)),'r');
   if fid == -1
   	disp(['READ_CSV_SLIT_FN.M: Error, ' filename(i,1:12) ' could not be opened.']);
   	x = -1; 
   	return
   end 
   
   % Preallocate memory for i'th row of meas_fn matrix:
   meas_fn(i,1:(2*wid + 1)) = zeros(1,(2*wid + 1));
   
   j = 0;
   %not_eof = 1;
   fgetl(fid);
   fgetl(fid);
   while ~feof(fid)	% loop over # of functions in one csv file
      % The following code retrives the centre coordinate, integration time,
   	% and mean CCD temperature from the csv file:
   	
      dt = fscanf(fid,'%44s',1);	% skip date and time fields
      disp(['  '])
      disp(strcat('date & time fields:',dt))	% debugging
      for m = 1:1:7
   	   fscanf(fid,'%5f',1);	% skip number
   	   fscanf(fid,'%1s',1);	% skip comma
   	end
   	wv_cent(i) = fscanf(fid,'%5f',1);	% get wavelength centre
   	fscanf(fid,'%1s',1);						% skip comma
   	int_time(i) = fscanf(fid,'%5f',1);	% get integration time
   	fscanf(fid,'%1s',1);						% skip comma
   	fscanf(fid,'%5f',1);						% skip number
   	fscanf(fid,'%1s',1);						% skip comma
      mean_TCCD = fscanf(fid,'%5f',1);		% get mean CCD temperature
   
   	% Make wv_cent and cent globally available by defining new variables:
%   	global swc
%   	global scent
%   	swc = wv_cent;

		% Convert wv_cent to a pixel centre:
      temp = cent(i);
      cent(i) = polyval(wv_fn, wv_cent(i));
         
   	% and round off:
   	if mod(cent(i),1) < 0.5   
   	   cent(i) = cent(i) - mod(cent(i),1);
   	else
   	   cent(i) = cent(i) +1 - mod(cent(i),1);
		end
      
      continue = 1;
      % Check that centre is same as previous centre of same function:
      if (j > 0) & (cent(i) ~= temp) 
         continue = 0;
         disp(['cent not equal to previous cent'])
      end
   	% Validate temperature range:
      if (mean_TCCD < T_min) | (mean_TCCD > T_max) 
         disp(strcat('READ_CSV_SLIT_FNS2: mean_TCCD = ', num2str(mean_TCCD),...
            ', out of range!'));
         continue = 0; 
      end
   	% Pixel range validation:
   	%if (cent(i) < 1) | (cent(i) > n_pix) continue = 0; end

   	if continue		% If cent and mean_TCCD are in range
   	   % Now read in the i'th slit function. We will read in values in a window
   		% that is 2*wid+1 pixels wide, centred around the centre one (i.e. +/- wid 
   		% pixels of cent(i)).
      
   	   fgetl(fid);	% now we are on the first line of comma-seperated data values
         
         if (cent(i)-wid > 0)	%cent(i)-wid = the starting pixel of the slit fn
            y = mod((cent(i) - wid), 200);
            % debugging:
            %if i == 11
            %   disp(['y 	   = ' num2str(y)])
            %   disp(['cent(i) = ' num2str(cent(i))])
            %end
            z = (cent(i) - wid) - y;
            for n = 1:1:(z/200)	% skip lines
      			fgetl(fid);
   			end
            for n = 1:1:(y-1)
   		   	fscanf(fid,'%8f',1);		% skip number
   	   		fscanf(fid,'%1s',1);		% skip comma
   			end
            [t1 t2] = min([2*wid 201-y]);
            % debugging:
            %disp(['t1 = ' num2str(t1) ', t2 = ' num2str(t2)])
            
            if (z == n_pix-200) t2 = 0; end	% prevent lineshape from running off
            											% end of the ccd pixel array.
         else	
            % consider case where slit function runs off beginning of pixel array:
            y = 1;
            z = 0;
            t1 = cent(i) + wid;  
            t2 = 0;
         end
   		% Now read in lineshape values:   			
         temp_meas_fn = zeros(1,2*wid+1);
         for n = 1:1:t1
            temp = fscanf(fid,'%8f',1);
            temp_meas_fn(n) = temp + temp_meas_fn(n);
            %meas_fn(i,n) = temp + meas_fn(i,n);
   	   	fscanf(fid,'%1s',1);		% skip comma
   		end
   		if t2 == 2		
   	   	fgetl(fid);		% skip to next line
            if feof(fid) 
               % note: this should not happen; if it does, something isn't doing
               % what it's supposed to. Jul30
               display(...
               	'READ_CSV_SLIT_FNS: Oops, function has fallen off end of file!');
      		end
      		n = n + 1;
      		while n < (2*wid+1)
               temp = fscanf(fid,'%8f',1);
               
               % debuggin:
               %if i == 11
               %   disp(['n_fns = ' num2str(n_fns)])
               %   disp(['i     = ' num2str(i)])
               %   disp(['n     = ' num2str(n)])
               %   whos temp
               %end
               
               temp_meas_fn(n) = temp + temp_meas_fn(n);
               %meas_fn(i,n) = temp + meas_fn(i,n);
      			fscanf(fid,'%1s',1);		% skip comma
      		   n = n + 1;   
      		end
   		end	% 	if t2 == 2
         
         % Now check that spectrum is above the cutoff level and below the
         % saturation level:
         peak = max(temp_meas_fn);
         if (peak > counts_cutoff) & (peak < counts_sat)
            % Add this spectrum to the total so far:
            meas_fn(i,:) = meas_fn(i,:) + temp_meas_fn;
            % Bring the noise:
         	if slit_fn_op.noise
            	dc_noise = darkcurrent(mean_TCCD, cent(i) - wid,...
            	   cent(i) + wid);
            	bias_noise = bias(mean_TCCD, cent(i) - wid,...
            	   cent(i) + wid);
      			k = 1:1:(2*wid+1);
      			meas_fn(i,k) = meas_fn(i,k) - dc_noise(k) - bias_noise(k);
            	% cast thy negative values into the void
            	meas_fn(i,k) = max(meas_fn(i,k), 0);
         	end		% if slit_fn_op.noise
            j = j + 1;   % increment count of # of slit functions averaged from one csv file
         	sf_files.temp(i,j) = mean_TCCD;
         elseif (peak < counts_cutoff)
            disp(['	spectrum counts below cutoff, discarding...'])
         else
            disp(['	spectrum counts above saturation, discarding...'])
         end
            
         % Skip lines to get to the beginning of the next spectrum
         if t2 == 2
            disp('t2 == 2')
            skip_lines = 8 - (z/200);
            disp(['skip_lines =    ' num2str(skip_lines)])
            disp(['cent(i) =       ' num2str(cent(i))])
         else
            skip_lines = 10 - (z/200);
            disp(['skip_lines =    ' num2str(skip_lines)])
         	disp(['cent(i) =       ' num2str(cent(i))])   
         end
         for n = 1:1:skip_lines
      		fgetl(fid);
   		end
           
      else		% continue < 1, spectrum invalid for some reason or another		
         disp(['	discarding spectrum...'])
         % Skip ahead to the next spectrum:
         for n = 1:1:11
            fgetl(fid);
         end
      end	% if continue    

	end	% while ~feof(fid)

	% Do crude averaging of meas_fn:
   % (perhaps should do a weighted average, weight by area?)
   if j > 0
      
      meas_fn(i,:) = meas_fn(i,:)./j;
      sf_files.filename(i,:) = strcat(dir_path,filename(i,1:12));
   	sf_files.n_avg(i) = j; 	% no. of averaged functions from a single csv file
   	disp(['No. of spectra averaged to make slit function no. ' num2str(i) ...
            ' = ' num2str(j)])
   else
      disp(['	i = ' num2str(i)])	%debugging
   
      disp(['All spectra in file were invalid (' ...
            filename(i,1:12) ')'])
      disp(['		discarding file...'])
      n_fns = n_fns - 1;	% reduce total number of functions
      disp(['Total no. of functions now reduced to ' num2str(n_fns)])    
      sf_files.files_not_used(p,:) = filename(i,:);
      p = p + 1;
      filename(i,:) = [];		% remove filename from list
      i = i-1;
      if n_fns == 0
         disp('READ_CSV_SLIT_FNS2: Uh-oh, we''re all out of slit functions!')
       	disp('		Exiting program...')
        	x = -1;
        	return
      end
   end	% if j > 0

   fclose(fid); 
	i = i + 1;
end	% for i = 1:1:n_fns
   
i = 1:1:n_fns;
coord(i) = cent(i) - wid;	% coord(i) = starting pixel of each slit function
len(i) = 2.*wid + 1;			% len(i) = length of each slit function (these are
									% actually all the same).
x = n_fns;
