power of blackbody radiation
37 views (last 30 days)
Show older comments
I need to calculate power of blackbody radiation at given temperature for wavelength ranging from 0 to 20 micrometer. And also need to plot the power vs wavelength.
I have intensity as a function of wavelength (Lam)
I2 =(2*h*c)./((Lam.^5).*(exp((h.*c)./(k.*T.*Lam))-1))
I integrated, but the results shows NaN
Please help me
clear all;
clc;
c=3*10^8; % speed of light in vaccum
h=6.625*10.^-34; % Planck constant
k=1.38*10.^-23; % Boltzmann constant
T= 500; % Temperatures in Kelvin
Lam=(0.0:0.01:20).*1e-6;
I1 =(2*h*c)./((Lam.^5).*(exp((h.*c)./(k.*T.*Lam))-1))
plot(Lam,I1)
avg_I1 = I1(1:length(Lam)-1)+diff(I1)/2;
% Integrated results
A = sum(diff(Lam).*avg_I1)
1 Comment
Clinton Edwards
on 4 Oct 2020
Edited: Clinton Edwards
on 4 Oct 2020
There are a few errors and enhancements which would be considered.
1) The most critical error is that the the calculation of I1 is missing a c^2 term and only has c in the first parathetical term. Instead of (2*h*c) is should be I1 =(2*h*c^2).
2) Less critical is that the first four significant figures of "h" are 6.626 not 6.625.
My proposed corrections are implemented in the code below with some axis labels.
Hope this helps! Nice job.
clint
% Physical Constants
c=3*10^8; % Speed of Light
h=6.626*10.^-34; % Planck constant (m^2*kg/s)
k=1.38*10.^-23; % Boltzmann constant (m^2*kg)/(s^2*K)
Lam=(Lam1:0.01:Lam2).*1e-6; % meters
dLam = Lam(2) - Lam(1); % Delta Lambda
% Calculate Blackbody Values For Vector Lam
I1 =(2*h*c^2)./((Lam.^5).*(exp((h.*c)./(k.*T.*Lam))-1));
plotI1 = I1*1e-12;
plotLam = Lam/1e-6;
plot(plotLam,plotI1,'linewidth',4)
xlabel('Wavelength (um)')
ylabel('Spectral Radiance (KW/m^/sr/wavelength)')
grid 'on'
Answers (1)
David Goodmanson
on 23 Apr 2019
HI rabindra,
The expression has a problem at the first point, lambda = 0. But I1 is so tiny near the origin that you can safely drop lambda= 0 and start the array at .01.
Things will work better using 2*h*c^2 rather than the incorrect 2*h*c. Also, the integration is correct but you can do the same thing more quickly with A=trapz(Lam,I1)
0 Comments
See Also
Categories
Find more on Symbolic Math Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!