function x = resid(ss, sp, ref, cba_sp, cba_ref, abs_fit)
% function resid(ss, wv_sp, wv_ref, sp, ref)
%	For given stretch and shift (contained in vector ss), 
%	resid does the following (in this order):
%		- convert ss to new a,b,c
%		- remake wavelength grids
%		- rebin sp and ref onto the new grid
%		- divide sp by ref
%		- take negative log (natural log)
%		- subtract slowly varying poly component
%		- rebin the lab cross sections to the new wavelength grid
%		- do a linear least squares fit of -log(sp/ref) using the 
%			lab cross sections, retrieving slant columns
%			i.e. -log(sp/ref) = sum((x-sec) x (slant column))
%									(sum over absorber species)
%			The slant columns are the coefficients of the fit
%		- get the residuals of the fit
%		- calculate the standard deviation of the fit

%	LAST MODIFIED:	July 17 / 00
%
%	ARGUMENTS:
%		ss				ss(1) = shift, ss(2) = stretch
%		sp				spectrum
%		ref			reference spectrum
%		cba_sp		initial a,b,c values for sp, format: [c b a]
%		cba_ref		initial a,b,c values for ref, format: [c b a]
%		abs_fit		array giving the indices of absorber gases to include
%						in the slant column fit. E.g. to include all gases,
%							abs_fit = [1 2 3 4 5 6 7] (or abs_fit = [1:7])
%						Or, to exclude OClO but include everything else,
%							abs_fit = [1 2 4 5 6 7]
%
%	VARIABLES:
%		lab_cs		matrix for absorption cross sections binned onto wv_ref
%		sc				slant columns (vector)
%		n_abs_fit	length(abs_fit), i.e. the # of absorbers to include in the sc fit
%						Obviously, n_abs_fit <= n_abs
%
%	CALLS:
%		wave_cent	gets "centre of gravity" of a spectrum
%		ss_convert	changes a,b,c params to new values using shift & stretch values
%		wave_grid	creates wavelength grid using a,b,c parameters
%
%	OUTPUT:			The residuals of the fit (cal_sp - sp_mk)
%
%OOPS, bin doesn't work...not using it right now. (July 14)
%		bin				rebins a spectrum on to a new wavelength grid

%	COMMENTS:
% For now, ref will be shifted as well as sp. (in the Fortran code, there is a flag
% set in WAV_CAL that determines whether or not ref gets shifted. See WAV_CAL pg 4
% and SHIFT1 pg 7)
%
%	July 18	-	note that code doesn't yet handle case of absorber wavelength ranges
%					being less than the wavelength range of cal_sp.

global resid_plot_results	% temp

% Globals from INIT.M:
global n_abs n_pix abs_gas abs_ex data_ex

% Counter and record keeping for Marq method:
global ss_count ss_record

global deg_calib flag_ring
% wavlength (wv) and cross sections (sm) for absorber gases:
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

% Show user the current shift/stretch values:
disp(strcat('shift = ', num2str(ss(1)), ', stretch = ', num2str(ss(2))))

% Record shift/stretch values and increment global counter:
ss_count = ss_count + 1;
ss_record(ss_count,1) = ss(1);
ss_record(ss_count,2) = ss(2);

%n_pix = 2000;		% this is now set globally
interp_method = 'cubic';
n_abs_fit = length(abs_fit);

% Generate wavelength grids for sp and ref:
wv_sp = wave_grid(cba_sp, n_pix);			% make wavelength grid for sp
w_cent = wave_cent(sp, wv_sp);				% get weighted centre of sp
cba_sp_new = ss_convert(ss, cba_sp, w_cent);	% get new a,b,c (functions of shift & stretch)
wv_sp = wave_grid(cba_sp_new, n_pix);		% remake wavelength grid for sp

%wv_ref = wave_grid(cba_ref, n_pix);			% make wavelength grid for ref
%w_cent = wave_cent(ref, wv_ref);				% get weighted centre of ref
% For now, use the value of w_cent from sp (this is what the old Fortran
% program does):
cba_ref_new = ss_convert(ss, cba_ref, w_cent);	
wv_ref = wave_grid(cba_ref_new, n_pix);		% remake wavelength grid for ref

% *****************************************************************
% temp fix: cut out portion of wv_ref that is out of wv_o3's range:
%continue = 1;
%m = 0;
%while continue
%   m = m + 1;
%   if wv_ref(m) > 270.0370
%      disp(strcat('m = ', num2str(m)))
%      continue = 0;
%   end
%end
%z = 1:(m-1);
%wv_ref(z) = [];
%ref(z) = [];
% Note, this does not change values of ref in the calling program. This also
% alters the pixel/wavelength relation, although it is easy to recover by simply
% doing a shift of m pixels. 
% *****************************************************************

%eff_n_pix = length(wv_ref);

% Note: the following (ie using the wavelenth cutoffs) might be better off in 
% a little subroutine of its own:
% Apply wavelength cutoffs as set in INIT.M (in the global struct "data_ex"):
%whos wv_ref
m = n_pix; %eff_n_pix;
while wv_ref(m) > data_ex.high_cut
   wv_ref(m) = [];
   ref(m) = [];
	m = m - 1;
end
m = 1;
while wv_ref(m) < data_ex.low_cut
   m = m + 1;
end
clear temp_ref temp_wv_ref
temp_ref = ref;
temp_wv_ref = wv_ref;
eff_n_pix = length(wv_ref);
for j = m:1:eff_n_pix
   wv_ref(j - m + 1) = temp_wv_ref(j);
   ref(j - m + 1) = temp_ref(j);
end
wv_ref((eff_n_pix - m + 1):eff_n_pix) = [];
ref((eff_n_pix - m + 1):eff_n_pix) = [];
clear temp_ref temp_wv_ref
eff_n_pix = length(wv_ref);
if length(ref) ~= eff_n_pix
   disp('RESID.M: Error, coding error in applying wavelength limits to data.')
   keyboard		% give control to the keyboard
end

% Bandaid for interpolation falling of end of range:
while min(wv_sp) > min(wv_ref)
   [t1 t2] = min(wv_ref);
   wv_ref(1:t2) = [];
   ref(1:t2) = [];
end
while max(wv_sp) < max(wv_ref)
   [t1 t2] = max(wv_ref);
   wv_ref(t2:length(wv_ref)) = [];
   ref(t2:length(wv_ref)) = [];
end
% This ensures that the range of wv_sp exceeds the range of wv_ref. If
% the range of wv_ref is greater than that of wv_sp, the interpolating
% function (interp1) will produce NaN values at the ends of cal_sp, 
% resulting in an error when polyfit is called.

eff_n_pix = length(wv_ref);

% Rebin sp onto the wavelength grid of sp:
cal_sp = interp1(wv_sp, sp, wv_ref, interp_method);

%cal_sp = bin(wv_sp, sp, wv_ref, interp_method);

% Divide sp by ref:
cal_sp = cal_sp./ref;

% Take negative natural log:
for m = 1:1:length(cal_sp)
   if (cal_sp(m) <= 0) cal_sp(m) = 1.0; end
end
% In MatLab, taking the log of zero gives -Inf as the answer. This won't crash
% the program; the -Inf will simply propagate throughout, producing a -Inf in
% the final result. But we don't want this to happen, so we set cal_sp = 1 so
% as to get a log of zero. 
cal_sp = -log(cal_sp);

% Fit and subtract slowly varying polynomial (to model aerosol absorption):
p = polyfit(wv_ref, cal_sp, 3);
cal_sp = cal_sp - polyval(p, wv_ref);

% Up to this point: runtime = 0.33 seconds for n_pix = 2000
% NOTE: at this point, we have cal_sp ready to fit, and binned on wv_ref. But
% remember that wv_ref may not be n_pix elements long anymore - some may have
% been chopped off, shortening it.

% Rebin absorber cross sections onto wv_ref. It is assumed here that the range
% of the absorber wavelength grids exceeds that of wv_ref. Otherwise, there
% will be problems (no fix for them yet)
for j = 1:1:n_abs_fit
   switch abs_fit(j)
   case 1
      lab_cs(j,:) = interp1(wv_o3, sm_o3, wv_ref, interp_method);
   case 2
      lab_cs(j,:) = interp1(wv_no2, sm_no2, wv_ref, interp_method);
   case 3
      lab_cs(j,:) = interp1(wv_oclo, sm_oclo, wv_ref, interp_method);
   case 4
      lab_cs(j,:) = interp1(wv_o4, sm_o4, wv_ref, interp_method);
   case 5
      lab_cs(j,:) = interp1(wv_h2o, sm_h2o, wv_ref, interp_method);
   case 6
      lab_cs(j,:) = interp1(wv_bro, sm_bro, wv_ref, interp_method);
   case 7
      if flag_ring > 0
   		lab_cs(j,:) = interp1(wv_ring, sm_ring, wv_ref, interp_method);
		else
   		lab_cs(j,:) = interp1(wv_no3, sm_no3, wv_ref, interp_method);
      end
   end	% "switch abs_fit(j)"
end		% "for j = 1:1:length(abs_fit)"
% Note that the row indices of matrix lab_cs no longer necessarily 
% correspond to the usual indices of the absorber gases.

% If we wish to exclude parts of the cross section for any absorber gas from the
% global scd fit, then we set those elements of the lab_cs matrix equal to zero now. 
% The excluded regions are stored in the global struct variable "abs_ex" (see 
% INIT.M for the structure).
for j = 1:1:n_abs_fit
   m = 1;
   while wv_ref(m) < abs_ex(abs_fit(j)).low_cut
      lab_cs(j,m) = 0;
    	m = m + 1;   
   end
   m = eff_n_pix;	% use eff_n_pix, instead of n_pix, in case we have shortened the data
   while wv_ref(m) > abs_ex(abs_fit(j)).high_cut
      lab_cs(j,m) = 0;
      m = m - 1;
   end
   %figure(j + 10)
   %plot(wv_ref, lab_cs(j,:));
end

%length(wv_ref)
%wv_ref(1)
%wv_ref(length(wv_ref))
%wv_o3(1:5)

% Do least squares fit to retrieve slant columns:
M = zeros(n_abs_fit);
b = zeros(1,n_abs_fit);
for j = 1:1:n_abs_fit
   b(j) = sum(cal_sp.*lab_cs(j,:));
   zero_flag = 0;
   for k = 1:1:n_abs_fit
      M(j,k) = sum(lab_cs(j,:).*lab_cs(k,:));
      if M(j,k) == 0
         zero_flag = zero_flag + 1;
      end
   end
   if zero_flag == n_abs_fit
      disp('RESID.M: Uh-oh...one row of M was all zeros! Exit from script.')
      return
      %n_abs_fit = n_abs_fit - 1;
      %abs_fit(j) = [];	% note that this does not change the value of abs_fit
      						% in the script that calls RESID.M
   %   M(j,:) = [];
   end
end

%M
%b'

M
% debugging fix: normalize the matrix
norm = min(min(M));
M = M./norm;
disp('After normalizing:')
M
b'

sc = M\b'

%M = zeros(n_abs);
%b = zeros(1,n_abs);
%for j = 1:1:n_abs
%   b(j) = sum(cal_sp.*lab_cs(j,:));
%   for k = 1:1:n_abs
%      M(j,k) = sum(lab_cs(j,:).*lab_cs(k,:));
%   end
%end
%sc = M\b;

% Remake spectrum using retrieved slant columns and absorber cross sections:
i = 1:1:length(cal_sp);
sp_mk(i) = 0;
for j = 1:1:n_abs_fit
   sp_mk(i) = sp_mk(i) + sc(j).*lab_cs(j,i);
end

% Temporary: let user know which gases have been used in the fit:
disp('Gases used in SC fit by RESID.M:')
disp(abs_gas(abs_fit))

%M
%sc
%b
%sp_mk(1:6)
% Get residuals:
res = cal_sp - sp_mk;

% Calculate standard deviation of the fit:
temp_sum = sum((cal_sp - sp_mk).^2);
%std_dev = 

disp(strcat('Chi squared = ', num2str(temp_sum)))

%x = res;		% Return residuals as the function output
x = temp_sum;		% return chi^2 as function output

% Do some plots, see what we got:
if resid_plot_results
disp('plotting...')
figure(1);
clf
hold on
plot(wv_ref, cal_sp, 'b');
plot(wv_ref, sp_mk, 'g');
plot(wv_ref, res, 'r');
title('cal\_sp = blue, sp\_mk = green, residual = red');
figure(2)
plot(wv_ref, cal_sp, 'b');
title('cal\_sp = blue');
end % "if plot_something"