MATLAB: Matlab fminsearch curve fitting HELP!!

blackbodycurve fittingCurve Fitting ToolboxfitfminsearchMATLAB

I'm very new to Matlab and working on a code based on planck's body radiation plots. I'm trying to use fminsearch to estimate temperature from a set of three wavelength(x) and irradiance(y) data points. Temperature is the only parameter to be estimated, others are either xy data or constants. My code is as follows:
    c=2.997*10.^8; % m/s (speed of light)
h=6.6261*10.^-34; % J.s (planck's constant)
k=1.38*10.^-23;   % T/K (boltzmann's constant)
xdata=[5.006e-07 5.104e-07 5.192e-07];
ydata=[3.223e+12 3.424e+12 3.606e+12]; 
figure(1)
clf
h1=plot(xdata,ydata,'ko','MarkerFaceColor','k');
set(gca,'YLim',[0,0.5*10^13]);
set(gca,'XLim',[0.3*10^-6,0.7*10^-6]);
set(gca,'XTick',xdata)
xlabel('Wavelength (m)');
ylabel('Energy Density');
t=3000;
pred = @(x,T) (c.^2)*h.*pi.*2)./((x.^5).*(exp((h.*c)./(k.*x.*T))-1));
hold on
h2=plot(xdata,pred,'r-');
err = @(x,T) sum((((c.^2)*h.*pi.*2)./((x.^5).*(exp((h.*c)./(k.*x.*T))-1))-ydata).^2);
T = fminsearch(@(x) err(xdata,x(1)),[3000]);
Please let me know where I'm wrong as I can't figure it out. Thanks.

Best Answer

First, ‘pred’ had mismatched parentheses, so it would be good for you to check it carefully to be certain you entered it correctly. I made one change it it that allowed it to run. Otherwise, you had some statements in the wrong order, and you were not using ‘pred’ correctly in your ‘err’ function. You are estimating ‘T’ not ‘x’, so you have to write ‘err’ to do that, and then optimise it with fminsearch.
It would have helped if you had described what you want to do and put comments in you code to describe its function. It took me a bit of time to figure out that you wanted to estimate ‘T’.
This works, and seems to be producing the correct result:
c=2.997*10.^8; % m/s (speed of light)
h=6.6261*10.^-34; % J.s (planck's constant)
k=1.38*10.^-23;   % T/K (boltzmann's constant)
xdata=[5.006e-07 5.104e-07 5.192e-07];
ydata=[3.223e+12 3.424e+12 3.606e+12]; 
figure(1)
clf
h1=plot(xdata,ydata,'ko','MarkerFaceColor','k');
set(gca,'YLim',[0,0.5*10^13]);
set(gca,'XLim',[0.3*10^-6,0.7*10^-6]);
set(gca,'XTick',xdata)
xlabel('Wavelength (m)');
ylabel('Energy Density');
T=3000;
pred = @(x,T) (c.^2)*h.*pi.*2./((x.^5).*(exp((h.*c)./(k.*x.*T))-1));
hold on
err = @(T) sum((pred(xdata,T)-ydata).^2);
T = fminsearch(err,[3000]);
h2=plot(xdata,pred(xdata,T),'r-');
Related Question