
function x = get_cs(f_ring, w_amf, ss)
% get_cs (w_amf, ss)
%   This function reads in laboratory absorption cross-sections. It can also
%   read in the wavelength-dependent AMFs (w_amf = 1 -> read them in, 
%   w_amf = 0 -> don't read them in) and shift n' stretch (i.e. rescale) the 
%   wavelength range for each individual gas species (ss = 1 -> rescale,
%   ss = 0 -> don't rescale).

% BUGS: 	- 	file read error when reading full array lengths, len_crs(i)
%			  	temp fix is to use len_crs - 1 when reading in cross sections (June 15)
%			- 	O3 cross section doesn't seem to look right (compare with general
%				appearance of cross sec. on pg. 84 of Deb Fish's thesis). All the 
%				others look okay (except Ring, I don't know what it supposed to look
%				like). (June 15)

% ARGUMENTS:
%				f_ring	whether to use NO3 or Ring cross section
%				w_amf		0 -> don't use wavelength-dependent AMF's; 1 -> use them
%				ss			0 -> don't shift/stretch wavelengths; 1 -> shift/stretch them
%
% INPUT FILES:
%		o3___crs.ext(1)	O3 cross section vs wavelength
%		no2__crs.ext(2)	NO2 cross section vs wavelength
%		oclo_crs.ext(3) 	OClO cross section vs wavelength
%		o4___crs.ext(4) 	O4 cross section vs wavelength
%		h2o__crs.ext(5)	H2O cross section vs wavelength
%		bro__crs.ext(6)	BrO cross section vs wavelength
%   	ring_crs.ext(7)	Ring effect cross section vs wavelength
%		     OR
%		no3__crs.ext(7)   NO3 cross section vs wavelength (use either ring OR NO3)
%				( NOTE: ext(i) is a string array given the 3-character extension of these
%				filenames. Different extensions denote varying array lengths for the
%				cross sections.)
%		ss.use				contains shift/stretch parameters for wavelength rescaling
%
% CALLS:			
%
% CALLED BY: 	input.m

% COMMENTS:
%	cs_lengths needs to be called first so that the array lengths, len_crs(i), for 
%	each gas species are known.

global n_abs		% # of absorbers

global ext			% cross section file extensions
global len_crs		% length of cross section arrays
global col			% unit vertical column for scaling the cross sections
% wavelength-dependent AMFs:
global wv_amf		% wavelength
global o3_amf		% AMF for O3
global no2_amf		% AMF for NO2
% wavlength (wv) and cross sections (sm):
global wv_o3
global wv_no2
global wv_oclo
global wv_o4
global wv_h2o
global wv_bro
global wv_ring
global wv_no3
global sm_o3
global sm_no2
global sm_oclo
global sm_o4
global sm_h2o
global sm_bro
global sm_ring
global sm_no3

dir_cs = '..\input_files\';

if (~(w_amf == 0 | w_amf == 1) | ~(ss == 0 | ss == 1))	% check input arguments
   disp('GET_CS.M: Error, invalid input arguments')
   x = 0;
	return
end

% Preallocate memory for cross section arrays. This should speed up the code:
wv_o3 	= zeros(1,len_crs(1));
sm_o3 	= zeros(1,len_crs(1));
wv_no2	= zeros(1,len_crs(2));
sm_no2	= zeros(1,len_crs(2));
wv_oclo	= zeros(1,len_crs(3));
sm_oclo	= zeros(1,len_crs(3));
wv_o4		= zeros(1,len_crs(4));
sm_o4		= zeros(1,len_crs(4));
wv_h2o	= zeros(1,len_crs(5));
sm_h2o	= zeros(1,len_crs(5));
wv_bro	= zeros(1,len_crs(6));
sm_bro	= zeros(1,len_crs(6));
if f_ring > 0	% use Ring instead of NO3
	wv_ring	= zeros(1,len_crs(7));
	sm_ring	= zeros(1,len_crs(7));
else   
   wv_no3	= zeros(1,len_crs(7));
	sm_no3	= zeros(1,len_crs(7));
end
   
N_ABS_LOC = 3;		% This is the same as n_abs. However, because we have to 
						% manually type in the filenames here, it is better to 
                  % manually type this in as well. It would be more elegant
                  % if this were done a different way (change it later...) (June 15)
file(1,1:8) = 'o3___crs';
file(2,1:8) = 'no2__crs';
file(3,1:8) = 'oclo_crs';
file(4,1:8) = 'o4___crs';
file(5,1:8) = 'h2o__crs';
file(6,1:8) = 'bro__crs';
if f_ring > 0
   file(7,1:8) = 'ring_crs';
else
   file(7,1:8) = 'no3__crs';
end
for i = 1:1:N_ABS_LOC
   file(i,9:12) = strcat('.',ext(i,1:3));
end

% for debugging, cut down array lengths:
for i = 1:1:N_ABS_LOC
   len_crs(i) = len_crs(i) - 1;		% BUG: 
   											% It seems that using the full array length,
                                    % len_crs(i), generates a file reading error.
                                    % I don't know which file the program is 
                                    % reading when the error occurs. (June 15)
end

% Read in cross sections. For the full arrays, this takes a while (about a minute
% or two).
for i = 1:1:N_ABS_LOC
   fid = fopen(strcat(dir_cs,file(i,1:12)),'r');
   if fid == -1 	% file was not opened
   	disp(strcat('GET_CS.M: Error opening ',strcat(dir_cs,file(i,1:12)),' - file not opened.'))
   	x = 0;
   	return
   end
   len_crs(i) = 1;
   j = 0;
   continue = 1;
   switch i
   	case 1  
         while continue 
            j = j + 1;
         	wv_o3(j) = fscanf(fid,'%f11',1);
            sm_o3(j) = fscanf(fid,'%f11',1);
            sm_o3(j) = sm_o3(j)*col(i);
            z = fgetl(fid);
            if feof(fid) 
               continue = 0;
               len_crs(i) = j;
            end
         end
      case 2
         while continue
            j = j + 1;
            wv_no2(j) = fscanf(fid,'%f11',1);
            sm_no2(j) = fscanf(fid,'%f11',1);
            sm_no2(j) = sm_no2(j)*col(i);
            z = fgetl(fid);
            if feof(fid) 
               continue = 0;
               len_crs(i) = j;
            end
         end
      case 3
         while continue
            j = j + 1;
            wv_oclo(j) = fscanf(fid,'%f11',1);
            sm_oclo(j) = fscanf(fid,'%f11',1);
            sm_oclo(j) = sm_oclo(j)*col(i);
            z = fgetl(fid);
            if feof(fid) 
               continue = 0;
               len_crs(i) = j;
            end
         end
      case 4
         for j = 1:1:len_crs(i)
            wv_o4(j) = fscanf(fid,'%f11',1);
            sm_o4(j) = fscanf(fid,'%f11',1);
            sm_o4(j) = sm_o4(j)*col(i);
            z = fgetl(fid);
         end
      case 5
         for j = 1:1:len_crs(i)
            wv_h2o(j) = fscanf(fid,'%f11',1);
            sm_h2o(j) = fscanf(fid,'%f11',1);
            sm_h2o(j) = sm_h2o(j)*col(i);
            z = fgetl(fid);
         end
      case 6
         for j = 1:1:len_crs(i)
            wv_bro(j) = fscanf(fid,'%f11',1);
            sm_bro(j) = fscanf(fid,'%f11',1);
            sm_bro(j) = sm_bro(j)*col(i);
            z = fgetl(fid);
         end
      case 7
         if f_ring > 0
            for j = 1:1:len_crs(i)
            	wv_ring(j) = fscanf(fid,'%f11',1);
            	sm_ring(j) = fscanf(fid,'%f11',1);
               sm_ring(j) = sm_ring(j)*col(i);
               z = fgetl(fid);
            end
         else
            for j = 1:1:len_crs(i)
            	wv_no3(j) = fscanf(fid,'%f11',1);
            	sm_no3(j) = fscanf(fid,'%f11',1);
               sm_no3(j) = sm_no3(j)*col(i);
               z = fgetl(fid);
            end
         end	% end of "if f_ring > 0"
   end	% end of "switch i"
   fclose(fid);
end	% end of "for" loop

% Need something here to handle wavelenth-dependent AMFs. I will worry about this 
% later, since it involves re-binning the AMF arrays, which requires writing
% a linear and/or polynomial interpolation routine. (June 12)
if w_amf == 1
   disp(strcat('GET_CS.M: Sorry...this code does not yet handle wavelength dependentAMFs. The management apologies for any inconvenience.'))
end

% If ss == 1, we shift and stretch the wavelengths for each species by the amounts
% given in the file SS.USE:
if ss == 1
   fid = fopen(strcat(dir_cs,'ss.use'),'r');
   if fid == -1 	% file was not opened
   	disp(strcat('GET_CS.M: Error opening SS.USE - file not opened.'))
   	x = 0;
   	return
   end
   for i = 1:1:N_ABS_LOC
      % Read in shift and stretch parameters for species i:
      z = fscanf(fid,'%1d',1);
      shift = fscanf(fid,'%10f',1);
      z = fscanf(fid,'%s4',1);
      stretch = fscanf(fid,'%10f',1);
      z = fgetl(fid);
      switch i
      	case 1
            for j = 1:1:len_crs(i)
               wv_o3(j) = wv_o3(j).*stretch + shift;
            end
         case 2
            for j = 1:1:len_crs(i)
               wv_no2(j) = wv_no2(j).*stretch + shift;
            end
         case 3
            for j = 1:1:len_crs(i)
               wv_oclo(j) = wv_oclo(j).*stretch + shift;
            end
   		case 4
            for j = 1:1:len_crs(i)
               wv_o4(j) = wv_o4(j).*stretch + shift;
            end
         case 5
            for j = 1:1:len_crs(i)
               wv_h2o(j) = wv_h2o(j).*stretch + shift;
            end
         case 6
            for j = 1:1:len_crs(i)
               wv_bro(j) = wv_bro(j).*stretch + shift;
            end
         case 7
            if f_ring > 0
               for j = 1:1:len_crs(i)
               	wv_ring(j) = wv_ring(j).*stretch + shift;
               end
            else
               for j = 1:1:len_crs(i)
               	wv_no3(j) = wv_no3(j).*stretch + shift;
               end 
            end	% end of "if f_ring > 0"
      end	% end of "switch i"
   end	% end "for i = 1:1:7"
	fclose(fid);   
end	% end "if ss == 1"
   
x = 1;
