How can I plot BB spectrum with logarithmic x axis and y 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);

Answer 1

[David Goodmanson]

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

[Marina Papathanasiou]

Thank you very much David! You are right but i wanted to plot it as you mentioned plugging in f = c/lambda to get a result that's a function of lambda. But if you change the temperature at your code the loglog plot is also distorted.

[David Goodmanson]

Hi Marina,

I'm not sure what you mean in your last sentence. If you change the temperature, the loglog plot changes. It should, since it's a function of T as well as lambda. I'm not sure what you mean by 'distorted'.