% spec_gen
%	This script takes the output from the slit_fn script and uses it
%	to generate a spectrometer spectrum from an input spectrum. The
%	input spectrum is binned as a function of wavelength, so we must 
%	convert it to pixels.
%
%	Runtime: about 90 seconds with 2000 pixels and a smoothing (slit)
%	function about 120 pixels wide for each pixel on the CCD.
%	(NOTE: that the slit functions are 120 pixels wide does not mean
%	that most of that range is significant - a large portion of the 
% 	(e.g.) 120 pixels may be zero or pretty close to it.)
%
%	   * ONLY RUN THIS SCRIPT AFTER RUNNING SLIT_FN FIRST! *

%	LAST UPDATE: July 7 / 00    by James Anstey
%	INPUT FILES:
%		file containing the raw spectrum
%
%	VARIABLES:
%	Input variables:
%		s_fn			matrix containing slit functions, ie the output of slit_fn
%		centre		tells the centre coordinate of each slit function in s_fn
%		n_pix			# of pixels on our CCD
%		max_len		max length of each slit function
%						the matrix s_fn is of dimension (n_pix)x(max_len)  (rows x columns)
%		spec_in		array containing our input spectrum
%		wv_in			array containing input wavelengths (ie x-axis for spec)
%							(units = nanometers)
%		a, b, c		params for pixel/wavelength conversion
%							pixel # = p, wavelength = L
%							p = a + b*L + c*L^2
%	
%	Output and other variables:
%		spec			array containing the rebinned input spectrum	
%		spec_out		array containing the simulated spectrometer spectrum. The
%							index numbers give the pixel numbers.
%
%	--> To summarize: the output consists of a spectrum, spec_out, that is
%	the convolution of an input spectrum and a total slit function for our
%	spectrometer. The slit function is generated by the script SLIT_FN.M
%

clear spec_in wv_in spec spec_out 
clear a b c j fid

tic
load slitfunctions		% load slitfunctions.mat
								% This file contains the results of running 
                        % SLIT_FN.M - the variables s_fn, n_pix,
                        % centre and max_len
toc                        
                     
%n_pix = 2000

% Initialize. Use values of a & b from A_B.REF (for lack of any better ideas):
dir_input = 'c:\james\analyse\input_files\';
fid = fopen(strcat(dir_input,'a_b.ref'),'r');
if fid == -1 	% file was not opened
   disp('SPEC_GEN.M: Error opening A_B.REF - file not opened.')
   return
end
fgetl(fid);
a = fscanf(fid,'%f20',1);
b = fscanf(fid,'%f20',1);    % a,b from absolute reference
fclose(fid);

c = 0;	% only linear relation between pixels and wavelengths for now.

% Read in spectrum:
fid = fopen(strcat(dir_input,'no2__crs.420','r'));
%fid = fopen('c:\james\analyse\input_files\no2__crs.420','r');
if fid == -1
   disp('oops, error opening file no2__crs.420');
   return
end
continue = 1;
j = 0;
while continue 
   j = j + 1;
 	wv_in(j) = fscanf(fid,'%f11',1);
   spec_in(j) = fscanf(fid,'%f11',1);
   fgetl(fid);
   if feof(fid) 
     	continue = 0;
      spec_in_length = j;
   end
end
fclose(fid);

figure(2);
plot(wv_in, spec_in, 'b');
title('no2 cross section read in by spec\_gen.m');


%wv_in(j) = 450;
%spec_in(j) = 10;

clear j
j = 1:spec_in_length;

% Convert wavelengths to pixel numbers (NOTE: these pixel #'s aren't
% necessarily integers - that happens a few lines down):
pix(j) = a + b*wv_in(j) + c*wv_in(j).^2;
% BUT this probably doesn't give a 1 to 1 correspondence of wavelength to 
% pixel number right off the bat. So we need to re-bin the spectrum so as
% to get a spectrum that is a continuous function of pixel number:
spec = zeros(1,n_pix);
spec = interp1(pix, spec_in, [1:n_pix], 'cubic');

spec_out = zeros(1,n_pix);
% Here is where we need the output from SLIT_FN: s_fn, centre, and max_len
clear m n
for m = 1:1:n_pix
   for n = 1:1:max_len
      if (m - centre + n >= 1) & (m - centre + n <= n_pix)
         spec_out(m) = spec_out(m) + spec(m - centre + n).*s_fn(m, n);
      end
   end
end

% Also try using MatLab's convolution function:
%spec_conv = conv(
% Actually, I don't think the convolution function is going to work. See
% lab notebook, page 39. (July 12)

%oops, wrong algorithm.
%for m = 1:1:n_pix
%   for n = 1:1:max_len
%      if (m + centre - n >= 1) & (m + centre - n <= n_pix)
%      	spec_out(m) = spec_out(m + centre - n) + spec(m).*s_fn(m, n);
%   	end      
%   end
%end

% Plot results:
figure(3);
clf	% clear figure
hold on;
plot(spec, 'r');
%plot(spec_out, 'b');
title('Input spectrum (red) vs pixel number');

toc

figure(4);
clf
plot(spec_out, 'b');
title('Output spectrum (blue)');

% Normalize and plot all spectra on same graph:
z_in 	= max(spec_in);
z 		= max(spec);
z_out	= max(spec_out);
norm_height = max([z_in z z_out]);
n_spec_in 		= spec_in.*(norm_height/z_in);
n_spec 			= spec.*(norm_height/z);
n_spec_out	 	= spec_out.*(norm_height/z_out);

figure(5)
clf
hold on
plot(pix, n_spec_in, 'r');
plot(n_spec, 'k');
plot(n_spec_out, 'b');
title('spec\_in = red, spec = black, spec\_out = blue');
