Help with code to fit a BlackBody to data to find temperature

Hello,

I am looking to fit a blackbody curve to some data I have and determine the temperature. I'm getting stuck as I am getting an error saying I have too many input arguments. I think this is about line 23? But I'm struggling to correct it.

Here's what I have so far:

%fit BlackBody to data 

%constants needed for BB equation 
c=2.997*10.^8; 
h=6.6261*10.^-34; 
k=1.38*10.^-23;

%wavelgnths and flux densities 
w = [0.000000445,0.000000477,0.00000163,7.625*10^-7,0.000000806,0.00000122,0.00000215,6.231*10^-7,0.000000658,3.543*10^-7,0.000000365,4.392*10^-7,0.000000231,3.465*10^-7,5.468*10^-7,0.000000291,0.000000212,0.000000551,9.134*10^-7];
f = [1.03919*10^-27,1.05624*10^-27,6.5907*10^-28,8.72017*10^-28,8.54111*10^-28,7.23846*10^-28,5.81202*10^-28,1.00213*10^-27,9.66438*10^-28,1.13496*10^-27,1.12742*10^-27,1.04994*10^-27,1.2677*10^-27,1.18296*10^-27,1.03105*10^-27,1.3626*10^-27,1.31748*10^-27,1.06062*10^-27,8.11377*10^-28];

%plot figure with data 
figure
plot(w, f, 'o', 'MarkerFaceColor', 'b', 'MarkerSize', 5);
xlabel('Wavelength (m)');
ylabel('Flux Density (W/m^2/Hz)');

%BlackBody function 
T = 30000; %starting point for Temperature 
BB = @(lam,T) 2.*h.*(c.^2)./((lam.^5).*1./(exp((h.*c)./(lam.*k.*T))-1));

%get difference between BB and data, and find minimum
diff = @(T) sum((BB(w,T)-f).^2);
T = fminsearch(diff,T);

%add fit to plot with red line
hold on 
plot(w,diff(w,T),'r-');

Answer 1

[Mathieu NOE]

hello

seems to me this was a minor issue in the way you code BB, do not use at the same time T as an array and a variable in the BB function

also do not shadow the diff function but use a different variable name (changed to delta)

but more important to me , I beleieve there must be an error somewhere in the equations (the constants looks ok to me at first glance) because there is no T value in the range 10^2 to 10^5 that allows the curves to match and the difference is huge (factor 10^40 at least) - see the second part of the code - a very simple test

also I believe the optimisation would result in better results if you would first fit a clean function between w and f instead of using scattered points with noise.

%fit BlackBody to data 
%constants needed for BB equation 
c=2.997*10.^8; 
h=6.6261*10.^-34; 
k=1.38*10.^-23;
%wavelgnths and flux densities 
w = [0.000000445,0.000000477,0.00000163,7.625*10^-7,0.000000806,0.00000122,0.00000215,6.231*10^-7,0.000000658,3.543*10^-7,0.000000365,4.392*10^-7,0.000000231,3.465*10^-7,5.468*10^-7,0.000000291,0.000000212,0.000000551,9.134*10^-7];
f = [1.03919*10^-27,1.05624*10^-27,6.5907*10^-28,8.72017*10^-28,8.54111*10^-28,7.23846*10^-28,5.81202*10^-28,1.00213*10^-27,9.66438*10^-28,1.13496*10^-27,1.12742*10^-27,1.04994*10^-27,1.2677*10^-27,1.18296*10^-27,1.03105*10^-27,1.3626*10^-27,1.31748*10^-27,1.06062*10^-27,8.11377*10^-28];
%plot figure with data 
figure
semilogy(w, f, 'o', 'MarkerFaceColor', 'b', 'MarkerSize', 5);
xlabel('Wavelength (m)');
ylabel('Flux Density (W/m^2/Hz)');
%BlackBody function 
T = 300; %starting point for Temperature 
BB = @(lam,T) 2.*h.*(c.^2)./((lam.^5).*1./(exp((h.*c)./(lam.*k.*T))-1));
%get difference between BB and data, and find minimum
delta = @(x) sum((BB(w,x)-f).^2);
Topt = fminsearch(delta,T)

 
Exiting: Maximum number of function evaluations has been exceeded
         - increase MaxFunEvals option.
         Current function value: 0.000000 

Topt = 8.3010e+20

%add fit to plot with red line
hold on 
semilogy(w,BB(w,Topt),'r-');

●

% test BB  for a couple of T values 
%plot figure with data 
figure
semilogy(w, f, 'o', 'MarkerFaceColor', 'b', 'MarkerSize', 5);
xlabel('Wavelength (m)');
ylabel('Flux Density (W/m^2/Hz)');
%BlackBody function 
T = logspace(2,5,100); %starting point for Temperature 
hold on 
for k = 1:numel(T)
semilogy(w,BB(w,T(k)),'r-');
end

[Mathieu NOE]

just saw this Fex submission, if that can help you :

Blackbody Radiance Calculations - File Exchange - MATLAB Central

[Caelan]

This is extremely helpful, thank you so much! I see where I went wrong with the code. And yes, it seems there's something wrong with my data - I'll go back to the source and find what I've done wrong. I did have to convert to flux density, so I'm assuming I've made an error there. Thank you for helping!

[Mathieu NOE]

as always , my pleasure !