MATLAB: How to plot BB spectrum with logarithmic x axis and y axis

black body loglog axis

Hello everyone! I am relative new in matlab and I am trying to plot a black body spectrum using Planck's equation but I want both x axis and y to be in logarithmic scale. But every time that I try to plot them both in logarithmic scale the graph's shape is distorted..I checked my units because I want to convert the meters to μm but still no luck..I post below the code. Thank you in advance everyone, I hope someone will help me!
c=3*10^14; % speed of light in vaccum [μm/sec]
h=6.625*10.^-34*10^12; %  Planck constant(m^2kg/s)*10^12 --> [microns]
k=1.38*10.^-23*10^12; %   Boltzmann constant [m^2 kg s^-2 K^-1]*10^12 -->[microns]
T=50; % Temperatures in Kelvin
b = 1.5;
Lam=(0.7:0.1:10000);
I1 = ((2*h*c)./(Lam.^3)).* (1./(exp((h*c)./(Lam*T*k))-1));
plot(Lam,I1);
set(gca, 'XScale', 'log');
%%set(gca, 'YScale', 'log');
figure
loglog(Lam,I1);

Best Answer

Hi Marina,
I believe the plots are basically correct, at least as far as they go. The shape does get changed quite a bit with the loglog plot, but if you zoom in you can see that the values all match that of the linear plot.
However, this plot is misleading. It is what one would get by using the expression for I(f,T) and plugging in f = c/lambda to get a result that's a function of lambda. In a plot with lambda as the independent variable, the implication is that you are plotting I(lambda,T). However, it does not quite work that way. The correct expression for I(lambda,T) is
I = (2*h*c^2./lambda.^5)./(exp(h*c./(lambda*k*T))-1)
This plot makes much more sense. For example follows the Wien displacement law, which the original plots do not.
The following code is a linear and loglog example in mks.
h = 6.626e-34;
k = 1.381e-23;
c = 299792458;
T = 50;
lambda = linspace(1e-5,1e-2,10000);
I = (2*h*c^2./lambda.^5)./(exp(h*c./(lambda*k*T))-1);
fig(1)
plot(lambda,I)
xlim([1e-5,1e-3])
xlabel(' wavelength, m')
ylabel('intensity, W / (m^2 steradian )')
fig(2)
loglog(lambda,I)
xlabel(' wavelength, m')
ylabel('intensity, W / (m^2 steradian )')
Related Question