% READ.M 
% This script reads in parameters and data needed to do the analysis.
% LAST UPDATE: July 18th
% Summary:
%		-	get names of spectra files, and store in "file_nms"
%		-	read in values from the parameter file ANALYSE.TXT
%		-	read A_B.REF to get filename of absolute reference spectra
%			and store in "std_name". Also read in initial a & b parameters.
%		-	read first and last pixel, and elevation angle, from the 
%			absolute reference spectra. Validate pixel range against
%			user-defined range read in from ANALYSE.TXT
%		-	read in absorber gas cross section arrays, and generate a
%			cross section for Rayleigh scattering
%		-	find out which data spectrum has the elevation angle closest
%			to zero, and store the index in "lowest"
%		-	get AMF's for each spectrum (read and interpolate)

% BUGS: 

% The matlab code to replace ANALYSE.FOR will be created in four "modules":
%
%	input data -> wavelength calibration -> SCD retrieval -> output results
%
% This particular script file contains my first attempt at the first task 
% (data input) (June 7)
%
%REQUIRED DIRECTORY STRUCTURE (July 17):
% 													CONTENTS:								
%		Analyse									Parent directory
%			|
%			|___Input_code					M-files (ie code) for data input tasks
%			|
%			|___Input_files				Input files (eg analyse.txt, etc)
%			|		|
%			|		|___Data_files			Files containing the spectra to be analysed
%			|
%			|___Mat_files					directory for storing useful .mat files
%			|
%			|___Calib_code					M-files for calibration of spectra
%
%
%
%
%
%*****************************************************************************
%
% ARGUMENTS:	none
%
% INPUT FILES (in order of reading): 
%		list.out				names of files with spectra to be analysed
%								header info is first 5 lines
%								then comes the filenames (one per line)
%
%		analyse.txt			various paramters 
%
%		a_b.ref				name of absolute reference spectrum (line 1)
%								initial a,b for abs. ref spectrum (line 2)
%								initial a,b for spectrum to fit (line 3)
%
%		wfn_col.dat			scaling factors for the lab cross sections
%
%		cross section files - see GET_CS.M for a list of input files
%
%		amf.lis				lists names of files containing AMFs for 
%								zenith sky and direct absorption measurements
%
% OUTPUT FILES:
%		cs_o3.mat			files containing the wavelength and spectra array
%		cs_no2.mat			for each of the absorber gas cross sections. E.g.,
%		cs_oclo.mat			cs_o3.mat contains the variables wv_o3 and sm_o3. 
%		cs_o4.mat			Stored in Mat_files directory.
%		cs_h2o.mat
%		cs_bro.mat
%		cs_ring.mat
%		cs_no3.mat
%
% VARIABLES:
%		op						struct with flags for various user-defined options:
%			cs_std_len		1 if absorber cross section lengths are predetermined,
%								0 if they are read (by checking for eof in the get_cs function)
%			get_cs			1 to read in and save cross sections, 0 not to
%			save_cs			1 to save the cross sections in .mat files, 0 not to
%
% CALLS:	cs_lengths		reads lengths of lab cross section arrays
%			get_cs_OLD		reads in the cross section data, older version of get_cs
%			get_cs			reads in the cross section data
%			binary_search	does a binary search (used to interpolate AMF data)
%
% CALLED BY: driving script for the whole analysis
%
% COMMENTS:
%
%	July 17	-	takes 494 seconds (8 min) to run when reading in full absorber
%					cross section arrays without predetermined array lengths. Takes
%					about 4 seconds if we skip reading in the cross sections.
%
%*****************************************************************************

% Globals set in the file INIT.M:
global abs_gas		% names of absorber gases
global n_abs		% # of absorber gases
global n_pix		% # of pixels on the horizontal axis of our CCD
global dir			% directory pathnames (struct)
global filenames	% various filenames (struct)

% Other global variables:
% (note: there are more globals declared at the end of this file)
global deg_calib
global flag_ring

global lowest
global a_ref b_ref a_sp b_sp

global ext
% BUG: Arghh...if ext is made global, then MatLab makes it of type double. WHY?!?
global len_crs
global col
% for wavelength-dependent AMFs:
global wv_amf
global o3_amf no2_amf
% wavlength (wv) and cross sections (sm):
%global wv_o3
%global sm_o3
global wv_o3 wv_no2 wv_oclo wv_o4 wv_h2o wv_bro wv_ring wv_no3
global sm_o3 sm_no2 sm_oclo sm_o4 sm_h2o sm_bro sm_ring sm_no3
% values read from ANALYSE.TXT:
global first_pixel last_pixel

% Set options in struct "op" 
op.cs_std_len = 0;
op.get_cs = 0;
op.save_cs = 0;
% Display options:

% SHOULD PROBABLY MAKE SOME OF THE FOLLOWING GLOBAL AND INITIALIZE THEM
% IN THE DRIVING SCRIPT: --> it is done (July 17)
% Initialization of important variables 
%n_abs = 7;		% # of absorbers (ie gas species & ring effect)
%n_pix = 1999;	% # of pixels


%dir_ip = '..\input_files\';	% Directory path to the input files
%dir_data = strcat(dir_ip,'data_files\');	% Directory path to the data files
%dir_cs_op = '..\Mat_files\';	% Directory path for storing cross section variables

%dir.ip = '..\input_files\';	% Directory path to the input files
%dir.data = strcat(dir.ip,'data_files\');	% Directory path to the data files
%dir.cs_op = '..\Mat_files\';	% Directory path for storing cross section variables
% THAT IS, THE ONES ABOVE SHOULD BE INITIALIZED BY THE DRIVING SCRIPT...?


debug = 0;		% toggles output of ANALYSE.TXT values (0 -> no o/p, 1 -> o/p)
an_yes = 1;		% 1 -> read ANALYSE.TXT values, 0 -> don't

% Read in filenames from LIST.OUT
% LIST.OUT contains the names of files that store the data to be analysed.
fid = fopen(strcat(dir.ip,'list.out'),'r');
if fid == -1 	% file was not opened
   disp('READ.M: Error opening LIST.OUT - file not opened.')
   return
end
% Count the lines in list.out so that we know when to stop reading filenames:
lines = 0;	
while feof(fid) == 0
   z = fgetl(fid);
   lines = lines + 1;
end
% Go back to beginning of file and skip the first five lines:
frewind(fid);
for i = 1:1:5
   z = fgetl(fid);
end
		% NOTE ON MATLAB:
		% I don't know why, but if you try to do the following loop as a "for"
		% loop, MatLab doesn't like it - you get this error message:
		%		??? Subscripted assignment dimension mismatch.
		% when trying to increment a counter as part of a character matrix in 
      % a "for" loop. Anyways, using a "while" loop fixes it.
      %
% Read in the names of spectra to be analysed:
i=1;
while i < (lines - 6)
   file_nms(i,1:8) = fscanf(fid,'%8s',1);
   file_nms(i,9) = '.';
   file_nms(i,10:12) = fscanf(fid,'%3s',1);
   z = fgetl(fid);
	i = i + 1;   
end
% Store the number of spectra for future use:
n_sp = lines - 7;
fclose(fid);

% Store file names in a global struct:
filenames.file_nms = file_nms;

if an_yes == 1 

% Read in values from ANALYSE.TXT
% ANALYSE.TXT lists the purposes of each of these variables
fid = fopen(strcat(dir.ip,'analyse.txt'),'r');
if fid == -1 	% file was not opened
   disp('READ.M: Error opening ANALYSE.TXT - file not opened.')
   return
end
% analyse.txt should be in the same format as it is when used by ANALYSE.FOR
n_recal = fscanf(fid,'%5d',1);  		z = fgetl(fid); % Skip the rest of the line (z = dummy variable)
n_cal = fscanf(fid,'%5d',1); 			z = fgetl(fid);
first_pixel = fscanf(fid,'%5d',1);
last_pixel = fscanf(fid,'%5d',1);	z = fgetl(fid);
lo_calib = fscanf(fid,'%5d',1);
hi_calib = fscanf(fid,'%5d',1);	  	z = fgetl(fid);
for i = 1:1:5
   lo_fit(i) = fscanf(fid,'%5d',1);
   hi_fit(i) = fscanf(fid,'%5d',1);
   z = fgetl(fid);
end
flag_shst = fscanf(fid,'%s',1);  	z = fgetl(fid);
flag_excl = fscanf(fid,'%5d',1); 	z = fgetl(fid);
flag_oclo = fscanf(fid,'%5d',1); 	z = fgetl(fid);
flag_bro = fscanf(fid,'%5d',1);  	z = fgetl(fid);
flag_ring = fscanf(fid,'%5d',1); 	z = fgetl(fid);
write_xs = fscanf(fid,'%5d',1);  	z = fgetl(fid);
diagnose = fscanf(fid,'%5d',1);  	z = fgetl(fid);
for i = 1:1:7
   shst_crs(i) = fscanf(fid,'%5d',1);
end
z = fgetl(fid);
deg_calib = fscanf(fid,'%5d',1);		z = fgetl(fid);
deg_anal = fscanf(fid,'%5d',1);		z = fgetl(fid);
i = 1;
ext = 'zoo';	% NOTE: I don't understand why this is, but we must 
					% initialize ext to some three (no more than three!) 
               % letter word. Otherwise, by declaring it to be global
               % we make it of type double, by default.
while i < 8
   ext(i,1:3) = fscanf(fid,'%3s',1);    z = fgetl(fid);
	i = i + 1;   
end
file_fmt = fscanf(fid,'%s',1);		z = fgetl(fid);
data_fmt = fscanf(fid,'%s',1);		z = fgetl(fid);
flag_pix = fscanf(fid,'%5d',1);		z = fgetl(fid);
z1 = fscanf(fid,'%s',1);	z = fgetl(fid);
z2 = fscanf(fid,'%s',1);	z = fgetl(fid);
z3 = fscanf(fid,'%s',1);	z = fgetl(fid);
z4 = fscanf(fid,'%s',1);	z = fgetl(fid);
res_nm = [z1;z2;z3;z4];
clear z1 z2 z3 z4;
diff = fscanf(fid,'%5f',1);			z = fgetl(fid);

fclose(fid);

% **************************************************************************
% Debugging stuff:
if debug == 1
   disp('READ.M: values read so far:')
   disp(strcat('	n_recal = ',num2str(n_recal)))
   disp(strcat('	n_cal	  = ',num2str(n_cal)))
   disp(strcat('	first_pixel = ',num2str(first_pixel)))
   disp(strcat('	last_pixel  = ',num2str(last_pixel)))
   disp(strcat('	lo_calib    = ',num2str(lo_calib)))
   disp(strcat('	hi_calib    = ',num2str(hi_calib)))
   for i = 1:1:5
      disp(strcat('	lo_fit(',num2str(i),') = ',num2str(lo_fit(i)),'  hi_fit(',num2str(i),') = ',num2str(hi_fit(i))))
	end
	disp(strcat('	flag_shst	= ',flag_shst))
	disp(strcat('	flag_excl	= ',num2str(flag_excl)))
	disp(strcat('	flag_oclo 	= ',num2str(flag_oclo)))
	disp(strcat('	flag_bro  	= ',num2str(flag_bro)))
	disp(strcat('	flag_ring 	= ',num2str(flag_ring)))
	disp(strcat('	write_xs 	= ',num2str(write_xs)))
	disp(strcat('	diagnose 	= ',num2str(diagnose)))
	for i = 1:1:7
	   disp(strcat('	shst_crs(',num2str(i),') = ',num2str(shst_crs(i))))
	end
	disp(strcat('	deg_calib = ',num2str(deg_calib)))
	disp(strcat('	deg_anal  = ',num2str(deg_anal)))
	for i = 1:1:7
	   disp(strcat('	ext(',num2str(i),') = ',ext(i,1:3)))
	end
	disp(strcat('	file_fmt = ',file_fmt))
	disp(strcat('	data_fmt = ',data_fmt))
	disp(strcat('	flag_pix = ',num2str(flag_pix)))
	for i = 1:1:4
	   disp(strcat('	res_nm(',num2str(i),') = ',res_nm(i,1:12)))
	end
	disp(strcat('	diff = ',num2str(diff)))
end   
% ***************************************************************************

end % end of "if an_yes == 1"

% Read A_B.REF to get parameters a & b, and the filename of the absolute reference:
fid = fopen(strcat(dir.ip,'a_b.ref'),'r');	% Open A_B.REF
if fid == -1 	% file was not opened
   disp('READ.M: Error opening A_B.REF - file not opened.')
   return
end
% Read in name of absolute reference spectrum, and values of a & b:
std_name = fscanf(fid,'%s',1);  z = fgetl(fid);
a_ref = fscanf(fid,'%f20',1);
b_ref = fscanf(fid,'%f20',1);  z = fgetl(fid);   % a,b for absolute reference
a_sp = fscanf(fid,'%f20',1);
b_sp = fscanf(fid,'%f20',1);  z = fgetl(fid);	 % a,b for spectrum to fit
fclose(fid);

% Store std_name in a globally accessible struct:
filenames.std_name = std_name;

% Now check that pixel range of absolute reference spectrum is correct. That is, 
% first_pixel and last_pixel (read from ANALYSE.TXT) should be the same values
% as n1 and n2, read from the reference spectrum file (std_name). First we
% open the absolute reference file and read off the lowest pixel, highest pixel,
% and elevation angle:
if (file_fmt == 'old') | (file_fmt == 'OLD')
   temp = strcat(std_name(1:9),'p',std_name(11:12));
   % By forcing there to be a "p" at a particular position in the filename,
   % we ensure that we are reading a file of the expected format.
   fid = fopen(strcat(dir.ip,temp),'r');
   if fid == -1 	% file was not opened
   	disp(strcat('READ.M: Error opening [',temp,'] - file not opened.'))
   	return
   end
elseif (file_fmt == 'new') | (file_fmt == 'NEW')
   temp = strcat(std_name(1:5),'p',std_name(7:12));
   % By forcing there to be a "p" at a particular position in the filename,
   % we ensure that we are reading a file of the expected format.
   fid = fopen(strcat(dir.ip,temp),'r');
   if fid == -1 	% file was not opened
   	disp(strcat('READ.M: Error opening [',temp,'] - file not opened.'))
   	return
   end
else
   disp(' READ.M: Error, file_fmt must be either "old" or "new"')
   return
end
% Now read in n1, n2 and elevation angle:
if (data_fmt == 'old') | (data_fmt == 'OLD')
   n1 = fscanf(fid,'%5d',1);	% first pixel
   n2 = fscanf(fid,'%5d',1);	% last pixel
   elev_ang_ar = fscanf(fid,'%10f',1);
elseif (data_fmt == 'new') | (data_fmt == 'NEW')
	n1 = fscanf(fid,'%5d',1);	% first pixel
   n2 = fscanf(fid,'%5d',1);	% last pixel
   fgetl(fid); fgetl(fid);
   elev_ang_ar = fscanf(fid,'%10f',1);
else
   disp(' READ.M: Error, data_fmt must be either "old" or "new"')
   return
end
fclose(fid);
% Now check n1 and n2 against first_pixel and last_pixel:
if (n1 ~= first_pixel) | (n2 ~= last_pixel)
   if n1 ~= first_pixel
      disp(strcat('READ.M: Error in [',std_name,'], n1 is not equal to first_pixel')); 
      disp(strcat(' n1 = ',num2str(n1),', first_pixel = ',num2str(first_pixel))); 
   end
   if n2 ~= last_pixel
      disp(strcat('READ.M: Error in [',std_name,'], n2 is not equal to last_pixel'));
      disp(strcat(' n2 = ',num2str(n2),', last_pixel = ',num2str(last_pixel))); 
   end
   return
end
% Tell the user what the absolute reference spectrum is:
disp(strcat(' File containing absolute reference spectrum:	',std_name))


% Read laboratory absorption cross sections for the various gases.
% First read in scale factors (unit vertical columns) for the cross sections:
fid = fopen(strcat(dir.ip,'wfn_col.dat'), 'r');
%fid = fopen('..\input_files\wfn_col.dat','r');
if fid == -1 	% file was not opened
   disp('READ.M: Error opening WFN_COL.DAT - file not opened.')
   return
end
i = 1;
while i <= n_abs
   abs_name(i,1:8) = fscanf(fid,'%8s',1);
   col(i) = fscanf(fid,'%f10',1);
   z = fgetl(fid);
   i = i + 1;
end
fclose(fid);

% Find lengths of cross section arrays:
if an_yes == 1
% The "if an_yes == 1" is here only because calling cs_lengths before 
% making ext a string causes an error. Debugging stuff.

if op.get_cs
	if op.cs_std_len
   	z = cs_lengths;
		if z == 0 
   		disp('READ.M: Error reading in lab cross section array lengths.');
			return   
		end
		% Now read in the cross sections, along with wavelength-dependent AMFs if desired:
		z = get_cs_OLD(flag_ring,0,0);
	else
      z = get_cs(flag_ring,0,0);	
   	% new function doesn't require predefined lengths for the absorber cross section 
   	% arrays (as the old Fortran program did), but instead reads in lengths from 
   	% the absorber cross section files (ie checks for end of file as it reads).
   end	% "if op.cs_std_len"
   % Save the cross section arrays, because it takes a long time to read all of them in:
   if op.save_cs
      save cs_o3 		wv_o3 	sm_o3
      save cs_no2 	wv_no2 	sm_no2
      save cs_oclo 	wv_oclo 	sm_oclo
      save cs_o4 		wv_o4 	sm_o4
      save cs_h2o 	wv_h2o 	sm_h2o
      save cs_bro 	wv_bro 	sm_bro
      if flag_ring > 0 
         save cs_ring 	wv_ring 	sm_ring
      else
         save cs_no3 	wv_no3 	sm_no3
      end
   end 	% "if op.save_cs"
end		% "if op.get_cs"

end	% "if an_yes", debugging stuff

% Generate a cross section for Raleigh scattering:
% (note that resolution is 0.5 nm)
for i = 1:1:801
   wv_ray(i) = 299.5 + 0.5.*i;
   sm_ray(i) = 7.564e-4/((wv_ray(i)/1000.0).^4);
end


% Now we will check the elevation angles of each data spectrum to find out which one
% is the closest to zero. This information will be used later when we calibrate the
% spectra.
% First read in elevation angle from each of the files to be analysed:
for i = 1:1:n_sp
   % Open the i'th data file:
   if (file_fmt == 'old') | (file_fmt == 'OLD')
      clear temp;
      temp = strcat(file_nms(i,1:9),'p',file_nms(i,11:12));
      fid = fopen(strcat(dir.data,temp),'r');
      if fid == -1 	% file was not opened
   		disp(strcat('READ.M: Error opening [',temp,'] - file not opened.'))
   		return
   	end
   else 	% that is, (file_fmt = 'new') | (file_fmt == 'NEW')
      clear temp;
      temp = strcat(file_nms(i,1:5),'p',file_nms(i,7:12));
      fid = fopen(strcat(dir.data,temp),'r');
   	if fid == -1 	% file was not opened
   		disp(strcat('READ.M: Error opening [',temp,'] - file not opened.'))
   		return
   	end   
   end
   % Read in elevation angle:
   if (data_fmt == 'old') | (data_fmt == 'OLD')
      for j = 1:1:3 elev_ang(i) = fscanf(fid,'%10f',1); end
   else	% that is, (data_fmt == 'new') | (data_fmt == 'NEW')
      fgetl(fid); fgetl(fid);
      elev_ang(i) = fscanf(fid,'%10f',1);
   end
   fclose(fid);  
   % Check to see if the i'th spectrum's elevation angle is the closest to zero so far:
   z = 90.0;
   lowest = 1;
   if abs(elev_ang(i) - (90.0 - diff)) < z
      lowest = i;
      z = abs(elev_ang(i) - (90.0 - diff));
   end
end 	% end of "for i = 1:1:n_sp"
% Now "lowest" contains the index of the spectrum with elevation angle closest to 0.


% Now we want to read in AMFs (not wavelength dependent ones).
% The AMFs we use depend on the nature of the spectra - zenith or direct absorption.
% Zenith spectra filenames have a "z" in the fifth position of their filenames, 
% whereas direct absorption spectra have something else. Dark current measurements
% have a "d" in the fifth position, and we don't want to deal with those (they should
% not be part of the list of input files, LIST.OUT).
% First, open the file AMF.LIS, which lists the names of files that contain the 
% AMFs:
fid = fopen(strcat(dir.ip,'amf.lis'),'r');
if fid == -1 	% file was not opened
   disp(strcat('READ.M: Error opening AMF.LIS - file not opened.'))
   return
end
name_zen = fscanf(fid,'%s',1);  fgetl(fid);
name_star = fscanf(fid,'%s',1);
fclose(fid);
if (file_nms(lowest,5) == 'z') | (file_nms(lowest,5) == 'Z')		% if zenith sky spectra
   fid = fopen(strcat(dir.ip,name_zen),'r');
   if fid == -1 	% file was not opened
   	disp(strcat('READ.M: Error opening [',name_zen,'] - file not opened.'))
      disp(' (this file contains AMFs for zenith sky spectra)')
      return
	end
   n_amf_data = fscanf(fid,'%5d',1);  fgetl(fid);
   for i = 1:1:n_amf_data
      zen_ang(i) = fscanf(fid,'%5f',1);
      amf_o3(i) = fscanf(fid,'%5f',1);
      amf_no2(i) = fscanf(fid,'%5f',1);
      fgetl(fid);
   end
   fclose(fid);
else		% if stellar spectra (ie direct absorption)
   fid = fopen(strcat(dir.ip,name_star),'r');
   if fid == -1 	% file was not opened
   	disp(strcat('READ.M: Error opening [',name_zen,'] - file not opened.'))
      disp(' (this file contains AMFs for zenith sky spectra)')
      return
   end
   n_amf_data = fscanf(fid,'%5d',1); fgetl(fid);
   for i = 1:1:n_amf_data
      zen_ang(i) = fscanf(fid,'%5f',1);
      amf_o3(i) = fscanf(fid,'%5f',1);
      amf_no2(i) = fscanf(fid,'%5f',1);
      fgetl(fid);
   end   
   fclose(fid);
end	% end of "if file_nms..."
% We wish to have an AMF for each spectrum, but the AMF data that we just read in 
% my not have an AMF value at the correct SZA (90 - elev_ang) for each of our
% data spectra. For those spectra that lie above the range of SZA used by our
% forward model, simply assign the AMF for the highest SZA available. For the 
% spectra within the available range of SZA, do a linear interpolation between
% SZA values to get an AMF for the particular SZA of each of our data spectra. 
for i = 1:1:n_sp	% loop over # of data files
   sza = 90.0 - elev_ang(i);
   if (sza > 97) 	% Our forward model only calcs values up to 97 deg...?
      amf(i,1) = amf_o3(n_amf_data);	% assign AMF with highest SZA
      amf(i,2) = amf_no2(n_amf_data);	% ditto
   else
      n_a = binary_search(sza, zen_ang, n_amf_data, 1);
      amf(i,1) = ((zen_ang(n_a)-sza)*amf_o3(n_a-1)+(sza-zen_ang(n_a-1))*amf_o3(n_a))/(zen_ang(n_a)-zen_ang(n_a-1));
		amf(i,2) = ((zen_ang(n_a)-sza)*amf_no2(n_a-1)+(sza-zen_ang(n_a-1))*amf_no2(n_a))/(zen_ang(n_a)-zen_ang(n_a-1));
   end
   for j = 3:1:n_abs
      amf(i,j) = amf(i,1);	% use O3 AMFs in place of other gases (for now)
   end
end 	% end of "for i = 1:1:n_sp
% And do the same for the absolute reference spectrum:
sza = 90.0 - elev_ang_ar;
if (sza > 97) 	
   amf_ar(1) = amf_o3(n_amf_data);	% assign AMF with highest SZA
   amf_ar(2) = amf_no2(n_amf_data);	% ditto
else
   n_a = binary_search(sza, zen_ang, n_amf_data, 1);
   amf_ar(1) = ((zen_ang(n_a)-sza)*amf_o3(n_a-1)+(sza-zen_ang(n_a-1))*amf_o3(n_a))/(zen_ang(n_a)-zen_ang(n_a-1));
	amf_ar(2) = ((zen_ang(n_a)-sza)*amf_no2(n_a-1)+(sza-zen_ang(n_a-1))*amf_no2(n_a))/(zen_ang(n_a)-zen_ang(n_a-1));
end
for j = 3:1:n_abs
   amf_ar(j) = amf_ar(1);	% use O3 AMFs in place of other gases (for now)
end
   
% MARKER - ALL VARIABLES INTRODUCED UP TO THIS POINT ARE NOTED IN VAR_INDEX.TXT,
% EXCEPT FOR GLOBAL DECLARATIONS (WHICH MAY OR MAY NOT BE NOTED).




