% Demos iterative fit to computer-generated blackbody spectrum
% to determine the color temperature (in K) and the emissivity.
% Simulates photon noise and random variation in the temperature
% Uses the fitblackbody.m function.
% T. C. O'Haver, May 2008
format compact
global emissivity

% Wavelength in nm
wavelength=[200:50:1600];
Pnoise=.01; % Photon noise
Tnoise=20; % Temperature noise
T=5000+Tnoise.*randn(size(wavelength)); % Temperature with random variation
emissivity=.1;
radiance = emissivity.*1.19111E+16*wavelength.^(-5)./(exp(14380000./(wavelength.*T))-1);
radiance = radiance+sqrt(radiance).*Pnoise.*randn(size(wavelength)); % Adds simulated photon noise

% Perform an iterative fit using the FMINSEARCH and fitblackbody functions
start=3000;    % generic initial guesses
options = optimset('TolX',0.1);  % Determines how close the model must fit the data
tic
Temperature=fminsearch(@(lambda)(fitblackbody(lambda,wavelength,radiance)),start);
toc

% Compute a model and plot it (blue line) along with 
% the original data (red points)
model=emissivity.*1.19111E+16*wavelength.^(-5)./(exp(14380000./(wavelength*Temperature))-1);
plot(wavelength,radiance,'r.',wavelength,model,'b')
xlabel( 'Wavelength, in nm' )
ylabel('Radiance, Watts/cm2/sr')
emissivity  
Temperature
