function tfit(absorbance)
%  Self-contained command-line demonstration function that compares the
%  TFit method (line 66) to the single-wavelength (line 57), simple
%  regression (line 58), and weighted regression methods (line 61). The
%  syntax is tfit(absorbance), where 'absorbance' is the true peak
%  absorbance of an absorber with a Lorentzian spectral profile of width
%  'width' (line 29), measured with a spectrometer with a Gaussian spectral
%  bandpass of width 'InstWidth' (line 30), fractional unabsorbed stray
%  light level of "straylight' (line 32), photon noise level of 'noise'
%  (line 31) and a random Izero shift of 'IzeroShift' (line 33). Plots the
%  spectral profiles and prints the measured absorbances of each method in
%  the command window. You can change the parameters in lines 28-33.
%
% Example: Just type tfit(absorbance) where absorbance = .001 - 100
% Tom O'Haver, August, 2006. Modified for Octave, Oct. 2012
format compact
format short g
global z c

% Arrays:
% y = true transmission spectrum, without noise or broadening
% InstFunction = instrument function (modeled as a Gaussian)
% yobsd = noisy instrumentally broadened spectrum
% f = frequency coordinate vector
% x = wavelength coordinate vector

% Initial values of the user-adjustable parameters:
ArrayLength=128;   % Number of points in signal
width=10     % FWHM of absorption peak 
InstWidth=20  % FWHM of broadening function (spectrometer slit width)
noise=0.01    % Random noise level when InstWidth = 1
straylight=.01 % May be a scalar or a vector of length ArrayLength (Slider adjustable)
IzeroShift=.01 % Random shifts in the 100% T intensity due to background absorption

% Define frequency and wavelength coordinate vectors
x=[1:ArrayLength]';
j=[-ArrayLength/2:(ArrayLength/2)-1]';
f=(ArrayLength/2)-abs(j);

% Calculate noisy instrumentally-broadened transmission profile, yobsd,
% by convoluting true transmission spectrum y with instrument function
% InstFunction and adding noise.  
% Note:  To model gaussian absorption, change 'lorentzian' to 'gaussian'
TrueSpectrum=lorentzian(x,(ArrayLength/2),width);
y=10 .^ (absorbance .* (-TrueSpectrum));
fy=fft(y);
InstFunction=gaussian(f,0,InstWidth);  % define Gaussian instrument function centered on zero
fa=fft(InstFunction);
fy1=fy.*fa;            % Convolve the transmisison profile with the instrument function (InstWidth)
yobsd=real(ifft(fy1));  % by multiplying Fourier transforms and inverse transforming the result.
yo=yobsd./sum(InstFunction);
yobsd=straylight+yo+((noise/InstWidth).*randn(size(yo))).*sqrt(yo);   % Add simulated photon noise
yobsd=yobsd.*(1-straylight);
yobsd=yobsd.*(1+IzeroShift.*randn); % Random shifts in Izero from sample to sample after instrument is zeroed

% Conventional methods
SingleWavelengthAbsorbance=-log10(yobsd(ArrayLength./2));
SimpleRegression=TrueSpectrum\(-log10(yobsd));
Background=ones(size(y));
weight=y;
WeightedRegression=([weight weight] .* [Background TrueSpectrum])\(-log10(yobsd) .* weight);

% Curve fitting method
start=10; % Because of the very large dynamic range of absorbance, two start values are
if SingleWavelengthAbsorbance<1,start=1;end  % used to prevent stalling on local optima.
lam=fminsearch(@(lambda)(fitM(lambda,yobsd,TrueSpectrum,InstFunction,straylight)),start);
results=[absorbance SingleWavelengthAbsorbance SimpleRegression WeightedRegression(2) lam];
disp(' ')
disp('          True A    SingleW    SimpleR      WeightR      TFit')
disp(results)

% (Optional) Plot spectral profiles (uncomment next 6 lines)
plot(x,real(yobsd),'r.',x,real(y),'g',x,real(c)*z,'b',x,gaussian(x,ArrayLength/2,InstWidth),'m:'); 
text(5,1.32,'Green = Reference spectrum      Dotted Magenta = Instrument function'); 
text(5,1.25,'                    Red = Observed T     Blue = Fit to observed T'); 
xlabel('Wavelength'); ylabel('Transmission');
title(['True absorbance = ' num2str(absorbance) '    Abs.Width = ' num2str(round(10*width)/10)  '   Inst.Width = ' num2str(round(10*InstWidth)/10) '    straylight= ' num2str(round(1000*mean(straylight))/10) '%' ]);
axis([0 ArrayLength 0 1.4]);

function err = fitM(lam,yobsd,Spectra,InstFun,StrayLight)
% Fitting function for broadened absorption of any number of components
% yobsd =  observed transmission spectrum (column vector)
% Sprecta = reference spectra for each component, one component/column
% InstFunction = Instrument function or slit function. (column vector)
% StrayLight = fractional stray light (scalar or column vector)
% Example: 
% options = optimset('TolX',0.000001);
% absorbance=fminsearch('fitM',1,options,[.0123 .0102 .0123 .0147],[.5 1 .5 .2],[1 .5 .0625 .5],.01)
% yobsd, Spectra, and InstFunction must have same number of rows (wavelengths)
%  T. C. O'Haver, August 2006
global z
global c
A = StrayLight + (10 .^ -(Spectra*lam'));
fy=fft(A);
fa=fft(InstFun);
fy1=fy.*fa;                
z=real(ifft(fy1))./sum(InstFun);   
c = z\yobsd;
q = z*c;
err = norm(q-yobsd);

function g = gaussian(x,pos,wid)
%  gaussian(x,pos,wid) = gaussian peak centered on pos, half-width=wid
%  x may be scalar, vector, or matrix, pos and wid both scalar
%  T. C. O'Haver, 1988
% Examples: gaussian([0 1 2],1,2) gives result [0.5000    1.0000    0.5000]
% plot(gaussian([1:100],50,20)) displays gaussian band centered at 50 with width 20.
g = exp(-((x-pos)./(0.6005615.*wid)) .^2);

function g = lorentzian(x,position,width)
% lorentzian(x,position,width) Lorentzian function.
% where x may be scalar, vector, or matrix
% position and width scalar
% T. C. O'Haver, 1988
% Example: lorentzian([1 2 3],2,2) gives result [0.5 1 0.5]
g=ones(size(x))./(1+((x-position)./(0.5.*width)).^2);
