power of blackbody radiation

37 views (last 30 days)
rabindra khadka
rabindra khadka on 22 Apr 2019
Edited: Clinton Edwards on 4 Oct 2020
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
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'

Sign in to comment.

Answers (1)

David Goodmanson
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)

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!