Plot curves with different colors inside the for loop and calculation of integral with trapz

I would like to plot blackbody spectral emissive power against wavelength in different temperatures. For that reason, I assigned a vector variable for wavelength and I used a for loop for the different temperatures. My code is the following:

clc
clear
% Constants
c0=2.998*10^8; % Speed of Light in vacuum
h=6.626*10^-34; % Planck's constant [J*s]
k=1.387*10^-23; % Boltzmann's constant [J/K]
sigma=5670*10^-8; % Stefan-Boltzmann constant [W/(m^2*K^4)]
lamda=0.1*10^-6:0.001*10^-6:25*10^-6; % Thermal Radiation 0.1-100μm (0.4-0.7μm visible light) [m]
n=1; % refractive index
T(1)=0; % Surface temperature of the sun
c=c0/n; % Speed of light not in vacuum
hta=1./lamda; % wavenumber
v=c./lamda; % frequency

for i=0:2:10
  if i<10
      T=i*500;
  else
      T=5762;
  end
  Ev=2*pi*h*v.^3*n^2./(c0^2*(exp(h*v/(k*T))-1)); 
% Planck's law  spectral emissive power function of frequency

    Elamda=2*pi*h*c0^2./(n^2*lamda.^5.*(exp(h*c0./(n*lamda*k*T))-1));% Planck's law  spectral emissive power funtion of wavelength

    Eb=trapz(Elamda) % Total Emissive power
    Eb1=n^2*sigma*T^4

    hold on
    set(gca,'XScale','log','YScale', 'log') % set log scaling instead of default linear
    axis([0.1 25 100 10^8]) %set axis limits
    xlabel('\bf Wavelength \lambda/μm')
    ylabel('\bf Blackbody Spectral Emissive Power, E_b_\lambda/(W/(m^2*μm))')
    title('\bf \fontsize{12} Blackbody Emissive Spectrum')

    plot(lamda*10^6,Elamda./10^6, '-.b') % if loglog is used instead of plot and hold on is also set then matlab plots in linear mode instead of logarithmic

end

First of all, I would like to have a different color for each curve. I tried to create a for loop and assign a string variable, the sort name of each color, which ultimately was set inside the plot command, but without success.

Finally, the *Eb *and *Eb1 *should give more or less the same results. Because Eb is the area under the curve of Elamda (computed with trapz) and the Eb1 is the integration of the Elamda (the integration of spectral emissive power, so the total emissive power). The problem occurs because of the step in lamda (the wavelength). Can you please explain how can I overcome this problem?

Thank you very much, Giorgos

Answer 1

[Matt Fig]

Use:

hold all

instead of:

hold on

And don't tell MATLAB to make every curve blue! You are telling MATLAB to make every curve blue with this string '-.b' when what you want is '-.' only.

As for your integrations, you can use the two argument form for TRAPZ:

trapz(x,y)

This will take care of the scaling problem. The real issue is that you forgot a decimal in sigma.

sigma=5.670*10^-8;

Putting this all together yields:

clc
clear
% Constants
c0=2.998*10^8; % Speed of Light in vacuum
h=6.626*10^-34; % Planck's constant [J*s]
k=1.387*10^-23; % Boltzmann's constant [J/K]
sigma=5.670*10^-8; % Stefan-Boltzmann constant [W/(m^2*K^4)]
% Thermal Radiation lamda 0.1-100?m (0.4-0.7?m visible light) [m]
lamda=0.1*10^-6:0.001*10^-6:25*10^-6; 
n=1; % refractive index
T(1)=0; % Surface temperature of the sun
c=c0/n; % Speed of light not in vacuum
hta=1./lamda; % wavenumber
v=c./lamda; % frequency
set(gca,'XScale','log','YScale', 'log') % set log scaling instead of default linear
axis([0.1 25 100 10^8]) %sett axis limits
hold all
xlabel('\bf Wavelength \lambda/?m')
ylabel('\bf Blackbody Spectral Emissive Power, E_b_\lambda/(W/(m^2*?m))')
title('\bf \fontsize{12} Blackbody Emissive Spectrum')
fprintf('\n\n%6s  |%15s%15s%10s\n','T','Eb','Eb1','%Diff')
fprintf([repmat('-',1,49),'\n'])

for ii=0:2:10
    if ii<10
        T=ii*500;
    else
        T=5762;
    end
    % Planck's law  spectral emissive power function of frequency
    Ev=2*pi*h*v.^3*n^2./(c0^2*(exp(h*v/(k*T))-1));
    % Planck's law  spectral emissive power funtion of wavelength
    Elamda=2*pi*h*c0^2./(n^2*lamda.^5.*(exp(h*c0./(n*lamda*k*T))-1));
    Eb=trapz(lamda,Elamda); % Total Emissive power
    Eb1=n^2*sigma*T^4;
    fprintf('%6i  |%15.2f%15.2f%10.2f\n',...
            T,Eb,Eb1,abs(Eb1-Eb)/max([Eb1,Eb,1])*100)
    plot(lamda*10^6,Elamda./10^6, '-.')
end
fprintf('\n\n')

[Giorgos Papakonstantinou]

Do you also know how can insert a legend with respect to the temperatures (in other words with respect to the T(i) inside the for loop)?

Thank you again for the help.

[Matt Fig]

legend('T=0','T=1000','T=2000','T=3000','T=4000','T=5762')

[Giorgos Papakonstantinou]

Thank you. I knew that already. I was hoping for a legend that would be affected automatically by the number of iterations.

It is peculiar that the result for Eb is 0 when T=0, because the temperature is in the denominator of the equation.

[Kennaa Gadaa]

pleas would you add legend

Answer 2

[Giorgos Papakonstantinou]

Thank you all for the replies and your remarks! The mistake about sigma was unacceptable... As for the trapz command, I didn't knew the two argument mode. You really saved me a lot of time.

Answer 3

[Sadanand Wachche]

Did anyone try to plot this in 3D? Taking T and lambda as on separate axis?