% PEAK_CAL.M
%	This script takes a spectrum and attempts to determine the pixel-wavelenth 
%	conversion, e.g. pixel # = p = a + b*wavelength. It is the first script
%	in a series of scripts used to calibrate the spectra.
% 		There will be an option for the user to choose a quadratic function as
%	well as the linear one, i.e. p = a + b*wavelength + c*wavelength^2. But this
%	is not in place now (June 23).
%		The calibration between wavelength and pixel number is determined by 
%	searching the spectrum for prominent spectral peaks of known wavelength, 
%  either four peaks from the Balmer series or some Fraunhofer lines for which
%	I don't know the source. 
%		This is the general algorithm: we take the spectrum and calculate the 
%	gradient spectrum. We then search for the points at which the gradient
%	cross zero (implying an extremum in the spectrum), thus giving us the peak
%  centres. We gradually narrow the search criteria until the desired number
%	of peaks is found, or until the routine fails and returns some error 
% 	message.

%  LAST UPDATE:	June 23

%	ARGUMENTS:
%
%	INPUT FILES:
%
%	CALLS:
%
%	CALLED BY:	driving script for the analysis program

% OTHER COMMENTS:
%	I will try to consistently use the subscript i to denote loops over pixel
% numbers. 
%	Local array names will be prefixed with the initials of the script, eg
% sp is locally copied into pc_sp (for peak_cal).
%	The first index of a matrix is always referred to as the row index, and the
% second index as the column index (which is, of course, the MatLab convention)
%	I will also try to remember to use lower case letters for vectors, upper
% case for matrices.

global first_pixel
global last_pixel

sp = [1 3 2 4 5 3 4 5 5 4 6 2 3 6 4 5 5 3 1 6 4 3 2 3 3 2 4 6 7 6 3 3 2];

n_reg = 1;	% TEMP

% POSSIBLY INITIALIZE THESE IN THE DRIVING SCRIPT INSTEAD OF HERE:
dir_ip = '..\input_files\';	% Directory path to the input files
dir_data = strcat(dir_ip,'data_files\');	% Directory path to the data files

% Peaks used for calibration:
calib_peak(1,1:4) =  [393.36   410.179  430.79   486.133];	% originally from w_lun
calib_peak(2,1:5) =  [379.79   383.539  388.905  397.007  410.179];	% from w_star 
calib_peak(2,6:10) = [434.047  486.133   0.0      0.0      0.0   ];	% from w_star

% Read in a,b for the reference spectrum. READ.M also reads these in. However,
% we do it again here just in case this script is called at a later point in 
% the program at which the file A_B.REF has been overwritten with new values
% of a & b.
fid = fopen(strcat(dir_ip,'a_b.ref'),'r');
if fid == -1 	% file was not opened
   disp('READ.M: Error opening A_B.REF - file not opened.')
   return
end
z = fgetl(fid);
pc_a_ref = fscanf(fid,'%f20',1);
pc_b_ref = fscanf(fid,'%f20',1);
fclose(fid);

% Initialize local array for the spectrum (THIS SHOULD PROBABLY BE DONE IN
% THE DRIVING SCRIPT):
nlo = 1;
nhi = length(sp);
%lo_n = first_pixel;
%hi_n = last_pixel;

% Get the log (natural log) of the spectrum:
if min(sp) <= 0
   disp('PEAK_CAL.M: Error, spectral data is zero or negative.');
   return
end
pc_sp = log(sp);
								% Do we _need_ to use the log spectrum if
								% all we are doing is searching for peaks?	

% Get gradient spectrum:
grad_sp = gradient(pc_sp);

% Degree of the fit. This should probably be read in from a file, ie set by the 
% user externally. The original program doesn't do that.
deg = 4;

% Adjust pixel ranges:	ACTUALLY...WHY DO WE NEED THIS? WE DON'T...I THINK.
%nlo = nlo + deg - 1;
%nhi = nhi - deg + 1;

% Fit a polynomial to the grad spectrum. This models the slowly varying part of
% the spectrum (e.g. aerosol scattering) so that we can subtract it off. We do
% a least squares fit; explanation is provided in Numerical Recipes (Chp 15).
% We formulate the matrix equation
%			A*b = y
% where b is a vector of the fit coefficients, of length deg+1. Y also has
% length deg+1, and A is of dimensions deg+1 rows by deg+1 columns.

% Initialize the matrix A:
%for j = 1:n_reg	% loop over regions, if we need to exclude regions of the
						% pixel range from the fit.
	% FUTURE CHANGE:
	%	For now, use MatLab sum fn to sum over pixels. If we need to exlude regions
	% of pixel space in the fit, can perhaps overload the sum fn to do this?
	% e.g.: call sum with an extra argument, n_reg. Check no. of args, and if
	% there is no extra arg, do normal sum. Otherwise, sum over each region and
	% then total the sums.

x = [nlo:nhi];	% create vector of the indices (ie pixel #s) of grad_sp

c = polyfit(x,grad_sp,deg);
		% Note that polyfit works just as well as the following matrix algebra.
      % However, I don't know if one would be able to control the sum function
      % within polyfit to accomodate different regions.
      
% **********************************************************      
% Comment out this code because polyfit does the same thing:
      
%for j = 1:(deg+1) 	% loop over rows
%   y(j,1) = sum(grad_sp.*(x.^(j-1)));	% equivalent to writing sum([argument],2)
%   for k = 1:(deg+1)	% loop over columns
%      A(j,k) = sum(x.^(j+deg-k));
%   end
%end

% Now solve for the coefficients:
%b = A\y;
% And change b into a row vector:
%b(1:deg+1,2) = 0;	% there's probably a more elegant way, but this works.
%b = b(1:deg+1);

% **********************************************************

% Debugging: display a plot to show what the polynomial and spectrum look like:
hold on		% so we can superimpose plots
plot(grad_sp,'r')
%plot(polyval(b,x),'b')
plot(polyval(c,x),'y')
hold off

% I will use, from now on, the polynomial fitted using polyfit.
% Subtract off the polynomial:
grad_sp = grad_sp - polyval(c,x);

% Now find the peaks in the spectrum:
cutoff = 6.0;
thresh = 6.0;
continue = 1;
tries = 0;		% logs # of attempts to find the peaks.

% stub for imitating spectra type selection. type = 1 -> w_lun, type = 2 -> w_star
type = 1;
n_peaks = 0;
if type == 1
   n_peaks = 4;
elseif type == 2;
   n_peaks = 10;
end

% TEMP STUFF:
continue = 0;
n_peaks = 1;

while continue == 1
   tries = tries + 1;
   % Calculate expected pixel centres of peaks:
   % two ways of doing same thing. leave them here now as example to myself.(June 23)
   j = 1:n_peaks; 
   exp_cent(j) = pc_a_ref + pc_b_ref*calib_peak(type,j);
   exp_cent
   % this one is cooler:
   exp_cent = pc_a_ref + pc_b_ref*calib_peak(1,:);
   exp_cent(:,n_peaks+1:length(exp_cent)) = [];
   exp_cent
   
   max_sp = max(grad_sp);
   min_sp = min(grad_sp);	% store these 'cause they're used several times
   
   peaks_so_far = 0;
   
   i = nlo;
   while (peaks_so_far < n_peaks) & (i < nhi)
      if (grad_sp(i) < min_sp/cutoff)
         width = 1;
         temp(width) = grad_sp(i);
         while width > 0
            % if we arrive at the right-hand side of a peak:
            if (grad_sp(i) >= max_sp/cutoff) & (grad_sp(i+1) < max_sp/cutoff)
               min_temp = min(temp);
               % if the minimum is low enough for this to be considered a peak
               			% HEY! why do we just look at the min? why not
               			% max as well? 
               if min_temp <= min_sp/thresh
                  peaks_so_far = peaks_so_far + 1;
                  % find the pixel # of the centre (stub)
                  cent(peaks_so_far) = (i - width)/2;
               end
               width = 0;
            else	% keep going to the right
               i = i + 1;
               width = width + 1;
               temp(width) = grad_sp(i);
               if (width > 20) width = 0; end
            end
         end	% end of "while width > 0"
      end		% end of "if (grad_sp(i) < min_sp/cutoff)"
      i = i + 1;
   end			% end of while "(peaks_so_far < n_peaks)"
   
   % Now check to see if we found the desired # of peaks:
%   if peaks_so_far < n_peaks
%      
%   end
   continue = 0;	

end	% end of "while continue == 1"

