Blackbody curve fitting to data

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

[Alex Sha]

Hi, Robbie, are the results below you want?

Root of Mean Square Error (RMSE): 1.77537844695845E-29
Sum of Squared Residual: 4.09755921890195E-57
Correlation Coef. (R): 0.959878177795913
R-Square: 0.921366116208802
Adjusted R-Square: 0.905639339450563

Parameter	Best Estimate
----------	-------------
emissivity	1.79392106216826E-23
temperature	41484.5441766574

[Image Analyst]

Only 13 points, and they don't even really fit the curve all that well? Can't you get any more points from your measuring device? Did you define wavelength??? What do you think reasonable values for emissitivy and temperature should be?

[Robbie Baillie]

Hi, thank you. The temperature should be around 25540K - but the fact you got an answer is a result - grateful for details on how you did it? Unfortunately I can't get any more results at the moment - I may be able to get a few more. They are the central wavelength's from different astronomy telescopes observing the AT2018COW explosion. I have actually converted the results from wavelength to frequency - i will amend this on the question. The idea is you fit the planck formula for different temperatures - the best fit gives you the temperature of the explosion. I then need to repeat this for another 50 odd time epochs to see how the Temperature changes over time.

[Amrtanshu Raj]

Hi,

Can you share the fitblackbody.m file? That would help us understand your problem better and propose a solution.

[Robbie Baillie]

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);

Answer 1

[Bjorn Gustavsson]

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