Blackbody curve fitting to data

17 views (last 30 days)
Robbie Baillie
Robbie Baillie on 20 Aug 2020
Answered: Bjorn Gustavsson on 26 Aug 2020
Hi, I am trying to fit a planck function to some data in order to determine the temperature. I have tried to adapt a previous programme done by O'Haver but I can't get either to work. I have tried the curve fitting app - writing the equation into custom equation - but that does not work either.
% Iterative fit to an experimentally measured blackbody spectrum
% to determine the color temperature (in K) and the emissivity.
% Uses the fitblackbody.m function.
% T. C. O'Haver, May 2008
format compact
global emissivity
% Enter data
Frequency=[6.97192E+14
6.46104E+14
1.83922E+14
3.71951E+14
2.41768E+14
1.36892E+14
4.55612E+14
5.44088E+14
3.33103E+14
8.21349E+14
1.33478E+15
1.15305E+15
1.55494E+15
]; % Frequency in Hz
radiance =[6E-29
5.689E-29
1.5439E-29
3.20102E-29
1.91032E-29
1.69931E-29
4.56719E-29
6.01172E-29
2.88265E-29
7.33041E-29
2.17793E-28
1.29134E-28
1.76482E-28
]; % Measured radiance in Watts/m-2Hz-1
% Perform an iterative fit using FMINSEARCH and fitblackbody.m
start=10000; % generic initial guess for blackbody temperature K
options = optimset('TolX',0.1); % Determines how close the model must fit the data
Temperature=FMINSEARCH('fitblackbody',start,options,frequency,radiance);
% Compute a model and plot it (blue line) along with
% the original data (red points)
model=emissivity.*1.474*E-50*frequency.^(3)./(exp(4.79924*E-11*frequency./(Temperature))-1);
plot(wavelength,radiance,'r.',wavelength,model,'b')
XLABEL( 'Wavelength, in nm' )
YLABEL('Radiance, Watts/cm2/sr')
emissivity
Temperature
  5 Comments
Amrtanshu Raj
Amrtanshu Raj on 26 Aug 2020
Hi,
Can you share the fitblackbody.m file? That would help us understand your problem better and propose a solution.
Robbie Baillie
Robbie Baillie on 26 Aug 2020
Hi, Thankyou for coming back to me.
function err = fitblackbody(~,frequency,y,~)
% Fitting function for a blackbody spectrum.
% T. C. O'Haver, May 2008
global emissivity
radiance=1.474E-50*frequency.^(3)./(exp(4.79924E-11./(frequency*Temperature))-1);
emissivity = radiance'\y';
z = radiance*emissivity;
err = norm(z-y);

Sign in to comment.

Answers (1)

Bjorn Gustavsson
Bjorn Gustavsson on 26 Aug 2020
This looks uggly. Matlab impresses me in this case, in that it still manages to obtain some kind of result. The uggliest problem here is that the error-function uses a global variable emissivity and assigns values to it. During the optimisation the values of this variable will change depending on what values of the search-parameters fminsearch tries. That cannot be a behaviour you want. I would make a modified version of fitblackbody:
function err = fitblackbody2(Temperature,frequency,radiance2fit,emissivity)
% Fitting function for a blackbody spectrum.
radiance = 1.474E-50*frequency.^(3)./(exp(4.79924E-11./(frequency*Temperature))-1);
% This scaling might be required to scale the radiances at different wavelengths
% radiance = radiance*emissivity;
err = norm(radiance2fit - radiance);
Then you'd have to modify the call to fminsearch:
Temperature = fminsearch(@(T) fitblackbody2(T,frequency,radiance,emissivity),start,options);
That way, you will have a fitblackbody2 function that doesn't rely on a global variable that changes values rather randomly (you cannot know what path fminsearch takes from the start-guess to the optimal parameter, so you should see the assignment of emissivity as an assignment to a random number).
HTH

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!