CalculateTfit(x,y,InstFunction,width,absorbance,InstWidth,straylight,noise,NoiseArray,IzeroShift)

function CalculateTfit(x,y,InstFunction,width,absorbance,InstWidth,straylight,noise,NoiseArray,IzeroShift)
% Performs calculation and plots graph for DemoTFit
% example:
% CalculateTfit(x,y,InstFunction,width,absorbance,InstWidth,straylight,noise,NoiseArray,IzeroShift)
global z
global c
% Define frequency and wavelength coordinate vectors
ArrayLength=length(x);
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);
% Unbroadened transmission profile, plotted as green line on the graph.
y=10 .^ (absorbance .* (-TrueSpectrum));
fy=fft(y);
% Broadening function , plotted as magenta dashed line
InstFunction=gaussian(f,0,InstWidth);  % defined Gaussian instrument function centered on zero
fa=fft(InstFunction);
fy1=fy.*fa;             % Convolution by multiplying Fourier transforms
yobsd=real(ifft(fy1));  % and inverse transforming the result.
yo=yobsd./sum(InstFunction);
NoiseArray=randn(size(yo));
% ybosd is the broadened transmission profile with noise, plotted as red dots on the graph.
% Use only one of the following two lines, to model either photon or constant noise
yobsd=straylight+yo+(noise.*NoiseArray).*sqrt(yo).*1/InstWidth;  % Adds simulated photon noise
% yobsd=straylight+yo+(noise.*NoiseArray).*1/InstWidth;  % Adds simulated constant noise
yobsd=yobsd.*(1-straylight);
yobsd=yobsd.*(1+IzeroShift.*randn); % Random shifts in Izero 
h=figure(1);

% Traditional methods
SingleWavelengthAbsorbance=-log10(yobsd(ArrayLength./2));  % Single-wavelength Beer's Law absornance
% Weighted regression 
Background=ones(size(y));
weight=y;
WeightedRegression=([weight weight] .* [Background TrueSpectrum])\(-log10(yobsd) .* weight);

% TFit method
options = optimset('TolX',0.000001);
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('fitM',start,options,yobsd,TrueSpectrum,InstFunction,straylight);
plot(x,real(yobsd),'r.',x,real(y),'g',x,real(c)*z,'b',x,gaussian(x,ArrayLength/2,InstWidth),'m:'); 
title([ 'Abs.Width = ' num2str(round(10*width)/10)  '   Inst.Width = ' num2str(round(10*InstWidth)/10) '    straylight= ' num2str(round(1000*mean(straylight))/10) '%   Noise = ' num2str(round(1000*noise)/10)  '%' ]);
text(1,1.07,'            Green = Unbroadened T      Magenta = Broadening function'); 
text(1,1.02,'                        Red = Observed T     Blue = Fit to observed T'); 
absorbancedisplay=round(10000*absorbance)/10000;
SingleWavelengthDisplay=round(10000*real(SingleWavelengthAbsorbance))/10000;
Regressiondisplay=round(10000*real(WeightedRegression(2)))/10000;
TFitdisplay=round(10000*lam(1))/10000;
xlabel([ 'Peak A = ' num2str(absorbancedisplay)  '   Single = ' num2str(SingleWavelengthDisplay) '    Reg. = ' num2str(Regressiondisplay) '  TFit = ' num2str(TFitdisplay)  ])