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

2012 9 29
3 个回答
回答已采纳
更新时间：2023 1 26
20 次查看（30 天）

0 个投票

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

请先登录，再回答此问题。

 采纳的回答

Matt Fig
Matt Fig 2012-9-29
编辑：Matt Fig 2012-9-29

0 个投票

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

4 个评论
显示 2更早的评论 隐藏 2更早的评论

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.
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.

请先登录，再进行评论。

更多回答（2 个）

0 个投票

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.

类别

在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by