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

for loopintegrationplot colortrapz

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

Best Answer

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