Fitting a blackbody curve with two variables - having issues with lsqcurvefit stopping and unsure how to use options effectively

Hello,

I am trying to fit a blackbody curve to some data I have and need to find the value for two variables. I have the following code:

format long g

%fit BlackBody to data

%constants needed for BB equation
% SI UNITS
c=2.998e8; 
h=6.6261e-34; 
kb=1.38e-23;


%wavelgnth (xdata) and flux densities (ydata). 

xdata = [0.000000212, 0.000000231, 0.000000291, 3.465E-07, 3.543E-07, 0.000000365, 4.392E-07, 0.000000445, 0.000000477, 5.468E-07, 0.000000551, 6.231E-07, 0.000000658, 7.625E-07, 0.000000806, 9.134E-07, 0.00000122, 0.00000163, 0.00000215];
ydata = [1.32933E-13, 1.11645E-13, 1.04652E-13, 6.34698E-14, 5.64255E-14, 5.39367E-14, 3.80465E-14, 3.66716E-14, 3.55176E-14, 2.93087E-14, 3.10419E-14, 2.40891E-14, 2.09837E-14, 1.42909E-14, 1.2889E-14, 1.01056E-14, 5.81715E-15, 3.5086E-15, 1.99137E-15];

%plot figure with data 
figure
plot(xdata, ydata, 'o', 'MarkerFaceColor', 'b', 'MarkerSize', 5);
xlabel('Wavelength (m)');
ylabel('Flux Density (wavelength depedence, W/m^3)');

%BB function - using variable names from MATLAB documentation for ease of
%learning to use fitting function. Here x(1) is the aplha multiplication
%factor and x(2) is temperature. xdata is wavelength

fun = @(x,xdata) x(1)*(2.*pi.*h.*(c^2)./(xdata.^5)).*(1./(exp((h.*c)./(xdata.*kb.*x(2)))-1));

%fitting

x0 = [1e-23 28000]; %starting point for fitting 

x = lsqcurvefit(fun,x0,xdata,ydata)

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

x = 1×2

1.0e+00 * 1.87324626094452e-30 28000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%plot fit 
hold on 
plot(xdata,fun(x,xdata),'r-');

The values produced for x show a value I would expect for x(2) but not for x(1). I would expect it to be around e-23. The solver is stopping, I have tried to edit the options for step count, but I'm unsure of what values to change within options and what would be reasonable numbers to use.

Anyone have any idea how to get a more accurate x(1) value and how to get the solver to run completely?

Thank you!

[Sam Chak]

Hi @Caelan

What estimation formula did you use to determine that the value for the parameter x(1) is expected to be on the order of ? Are all blackbody radiations characterized in this manner?

Answer 1

[John D'Errico]

Your data does not seem to fit terribly well to the model you propose.

c=2.998e8; 
h=6.6261e-34; 
kb=1.38e-23;
%wavelgnth (xdata) and flux densities (ydata). 
xdata = [0.000000212, 0.000000231, 0.000000291, 3.465E-07, 3.543E-07, 0.000000365, 4.392E-07, 0.000000445, 0.000000477, 5.468E-07, 0.000000551, 6.231E-07, 0.000000658, 7.625E-07, 0.000000806, 9.134E-07, 0.00000122, 0.00000163, 0.00000215];
ydata = [1.32933E-13, 1.11645E-13, 1.04652E-13, 6.34698E-14, 5.64255E-14, 5.39367E-14, 3.80465E-14, 3.66716E-14, 3.55176E-14, 2.93087E-14, 3.10419E-14, 2.40891E-14, 2.09837E-14, 1.42909E-14, 1.2889E-14, 1.01056E-14, 5.81715E-15, 3.5086E-15, 1.99137E-15];

I'll use my own fminspleas to do the fitting, as it is far more robust to problems.

fun = @(x,xdata) (2.*pi.*h.*(c^2)./(xdata.^5)).*(1./(exp((h.*c)./(xdata.*kb.*x))-1));

[INLP,ILP] = fminspleas({fun},28000,xdata,ydata)
INLP =
          14620.5497503281
ILP =
      1.48319389723619e-29
      
plot(xdata,ydata,'o')
hold on
plot(xdata,ILP*fun(INLP,xdata),'-s')

Sigh. The plot did not come out terribly clearly. But you can see it chose a black body temperature around half of what you expected.

[Caelan]

This was really helpful thanks! You're right, I had an error in the data, I have rectified this now and the code worked perfectly.

[Caelan]

Also - since this is your own function - can I ask if there is a way to determine the errors on the fitted INLP and ILP values?