# Planck Laws

### Raveau david (view profile)

06 Oct 2005 (Updated )

We cal the number of photon of a black bodies.

planck4.m
%    Calculate the number of photon with two methods in a scale of wavelength and
%    the curves for a blackbody

%%   Define the wavelength
%%   Range and increment.

%We calcul the number of photon with two methods.
%We dertermine the differnce of luminance for two (or many ) wavelength.
%We can find the error with the temperature and the wavelength
%for different acquisition time and many geometry.

% David Raveau, Phd
% david.raveau-etudiant@ineris.fr
% Setember 2005.

clear all;
close all;
t0 =cputime;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Ce programme permet de calculer le nombre de photon et la luminance pour une longueur d onde suivant diffrentes temprature')

%Choix de la determination de la luminance en fonction de lambda et de la temprature choix 1 methode intgrale ou choix 2 par courbe rduite de planck
%%%%%%%%  Define a few constants for two unity
c = 2.997930 * 10^8;             %   Speed of light
h = 6.6245*10^(-34);     %   Planck's constant
k = 1.38033*10^(-23);      %   Boltzmann's constant
c1 = 2*h*(c*1e6^2)^2;    %   First constant in Planck formula en W/m2/m4/sr
c2 = h*c/k *1e6;     %  Second constant in Planck formula en m.K
C1 = 2* h* c^2; %First constant Planck W/m2/sr
C2 = h*c/k; % Second constant Planck m.K
b = (c1  / (2898^5 *(exp (c2 / 2898)-1))); %constant b for the calcul of small curvepour
Sigma = b * c2 / (c2/2898.8) * 1.5203 ; % Constatn of Stefan or  K = 5.67e-8 i

%First iteration for many calcul
boucle = 0;
Nombre = 1;
l=0; %initialisation de la boucle
donnees = {'Valeur'};
titre = 'Nombre de calcul sur une Longueur d onde';
defaut = {num2str(Nombre)};
nom = inputdlg(donnees,titre,1,defaut);
if size(nom) ~= 0
Nombre = str2num(nom{1});
end
Resultat = 0;
n = 0 ; % Initialisation de la boucle

for boucle = 1: Nombre;
% Value type
Tmin = 773; % Temperature min
Tmax = 1673; % Temperature max
DeltaT = 20; % incrementation of temperature
Lambda = 0.5; % Choise of wavelength center in m
DeltaLambda = 0.01; % Choise of value of incrementation for the wavelength
Lambdamin = 0.3; %wavelength min
Lambdamax = 5; % wavelength max
Tau = 1; %It is the maxima of transmision on the wavelength center
%windows for temperature and wavelength and tau
donnees = {'Temprature min en K','Temprature max en K','Delta Temprature en K','Longueur d onde en m ',...
'Delta Longueur d onde en m ','Longueur d onde mini en m ','Longueur d onde maxi en m ','transmission  Lambda'};
titre = 'Donnes pour le calcul de la luminance entre deltalambda en fonction de la temprature :';
defaut = {num2str(Tmin),num2str(Tmax),num2str(DeltaT),num2str(Lambda),num2str(DeltaLambda),...
num2str(Lambdamin),num2str(Lambdamax),num2str(Tau)};
nom = inputdlg(donnees, titre, 1, defaut);
if size(nom) ~= 0
Tmin = str2num(nom{1});
Tmax = str2num(nom{2});
DeltaT = str2num(nom{3});
Lambda = str2num(nom{4});
DeltaLambda = str2num(nom{5});
Lambdamin = str2num(nom{6});
Lambdamax = str2num(nom{7});
Tau = str2num(nom{8});
end
% Value type
Temps_d_acquisition = 1e-6; %seconde
Diametre_objet= 0.024 ; %metre
Distance_lentille_detecteur = 0.006; % metre
Distance_lentille_objet= 0.072;
Diametre_lentille = 0.006; %metre
Diametre_fibre = 0.2e-3;
%windows for dtermine the solid angular and the time for the acquisition
donnees2 = {'Temps d acquisition en s','Diamtre de la surface metrice en m','Distance entre objet et la lentille en m ',...
'Diamtre de la lentille en m'};
titre2 = 'Donnes pour le calcul de l angle solide, de la surface et du temps:';
defaut2 = {num2str(Temps_d_acquisition),num2str(Diametre_objet),num2str(Distance_lentille_objet),...
num2str(Diametre_lentille)};
nom2 = inputdlg(donnees2, titre2, 1, defaut2);
if size(nom2) ~= 0
Temps_d_acquisition = str2num(nom2{1});
Diametre_objet = str2num(nom2{2});
Distance_lentille_surface = str2num(nom2{3});
Diametre_lentille = str2num(nom2{4});
end
%Calul of aera and angular
Surface_vise = Distance_lentille_surface^2 / Distance_lentille_detecteur^2 *(pi * (Diametre_fibre^2)/4)
Diametre_vise = sqrt(Surface_vise/pi*4);
Surface_objet = pi * (Diametre_objet^2)/4

Angle_solide = pi * (Diametre_lentille / 2 ) ^ 2 / Distance_lentille_surface^2 % Stradiant

Surface_mesure = 0;

if Diametre_vise < Diametre_objet

Surface_mesure = Surface_vise

else

Surface_mesure = Surface_objet

end

% Dtermination des deux vecteur Temprature et Longueur d'onde et initialisation
T = 0;
T = Tmin: DeltaT :Tmax;
nbr = length(T);
Longueur_d_onde = 0 ;
Longueur_d_onde = linspace(Lambdamin, Lambdamax, nbr); % Le vecteur longueur d'onde a la mme taille que T
D_Longueur_d_onde = (Lambdamax - Lambdamin )/ nbr;

%Creation de la gausienne autour du point de Lambda  pour trouver les bords des deux ailes
Longueur_d_onde1 = 0 ;
Tau_lambda1 = 0;
I = 0;
Lambda_a = 0;
Lambda_b = 0;
Precision = 0.0001; %Valeur pour la dtermination du calcul de lambda_a et Lambda_b
Longueur_d_onde1 = linspace(Lambda - 2 * DeltaLambda , Lambda + 2 * DeltaLambda , nbr); % ce vecteur permet de tracer la gaussienne et de plus de trouver les deux points a et b de ce tte meme gaussienne
Tau_lambda1 = exp(-1/2*((Longueur_d_onde1 - Lambda) / ( DeltaLambda / 2.35)) .^2 ) * Tau;
I = find(Tau_lambda1 > Precision);
Lambda_a = Longueur_d_onde1 (I(1));
Lambda_b = Longueur_d_onde1 (I(end));

% Creation de la gaussienne mais cette fois on le recalcul sur une plage de longueur d'onde compris entre lambda a et lambda b%
Longueur_d_onde2=0;
Tau_lambda2 = 0;
Longueur_d_onde2 = linspace(Lambda_a , Lambda_b , nbr ); % c'est ce nouveau vecteur de logueur d'onde qui doit servir au diffrente intgration entre les deux longueurs d'onde
Tau_lambda2 = exp(-1/2*((Longueur_d_onde2 - Lambda) / ( DeltaLambda / 2.35)) .^2 ) * Tau;
D_Longueur_d_onde2 = (Lambda_b - Lambda_a )/ nbr;
Longueur_d_onde3 = 0;
Longueur_d_onde3 = linspace(Lambda - (DeltaLambda/2), (Lambda + (DeltaLambda/2)) , nbr ); % on a un  nouveau vecteur d'onde mais cette fois ci compris sur delta de longueur d'onnde

%on a donc un veteur Temprature de mme taille que le vecteur longueur d'onde ainsi que celui compris entre Lambda a Lambda b
%On calcule l'nergie photonique
Energie_Photon_Lambda2 = h * c ./ (Longueur_d_onde2 * 1e-6); %On se place en m pour les longueur d'onde car on souhaite avoir des Joules
Energie_Photon_integer = trapz(Longueur_d_onde2,Energie_Photon_Lambda2)*1 / ((Lambda_b-Lambda_a)); %On integre sur delta lambda pour avoir l'energie moyenne
Energie_Photon_integer2 = 1/((Lambda_b-Lambda_a)*1e-6) * h *c *log(Lambda_b/Lambda_a);%On integre sur delta lambda pour avoir l'energie moyenne

%Initialisation des vecteurs
l = 2*n+1;
Luminance = 0;
Luminance_DeltaLambda = 0;
Luminance_DeltaLambda_Tau = 0;
% la constante nous donne le l'angle solide  le temps de mesure et la gomtrie de la surface
cte = Angle_solide * Temps_d_acquisition * Surface_mesure;
% Creation de la boucle sur la temprature frquence ety lngueur d'onde.
for i =1:nbr

for j=1:nbr

Luminance(j,i) = c1*Longueur_d_onde(j)^(-5) / ( exp(c2/(Longueur_d_onde(j)*T(i))) -1) ; % exact loi de planck

Luminance_DeltaLambda(j,i) = c1*Longueur_d_onde2(j)^(-5) / ( exp(c2/(Longueur_d_onde2(j)*T(i))) -1) ;%on calcul la luminance exacte sur les deux borbs des ailes de la gaussiennes

Luminance_DeltaLambda_Tau(j,i) = Luminance_DeltaLambda(j,i) * Tau_lambda2(j);

Luminance_DeltaLambda3(j,i) = c1*Longueur_d_onde3(j)^(-5) / ( exp(c2/(Longueur_d_onde3(j)*T(i))) -1);%on calcul la luminance exacte mais sur Lambda +ou- deltalambda

end

Luminance_integer(i) = trapz(Longueur_d_onde2, Luminance_DeltaLambda(:,i)); % ceci est l'inetrgale compris entre Lambda a et lambda b sans avoir tenucompte de la gaussienne
Luminance_integer_tau(i) = trapz(Longueur_d_onde2, Luminance_DeltaLambda_Tau(:,i));
Resultat(i,l) = T(i) ;%dans cette matrice on garde en mmoire sur la premire colonne le vecteur T et sur la seconde
Resultat(i,l+1) = Luminance_integer(i); % en incrmentant la luminance sur deltalambda % si on souhaite avoir la luminance en fonction de Tau il suffit d'ajouter tau  la fin de la ligne

Nombre_Photon(i) = Luminance_integer_tau(i) /  Energie_Photon_integer *cte ; % Nombre de photon s-1 m-2 et sr-1 c'est pour cela
%
Luminance_integer3(i) = trapz(Longueur_d_onde3, Luminance_DeltaLambda3(:,i));

Luminance_Lambda(i) = c1*(Lambda^(-5)) / ( exp(c2/(Lambda*T(i))) -1)*(DeltaLambda) ; % exact loi de planck  Lambda

Nombre_Photon_Lambda(i) = Luminance_Lambda(i) / (h*c/(Lambda*1e-6))*cte; % ce calcul permet une premire aproximation

end
%Representation des diffrents graphiques
figure
plot(Longueur_d_onde2 , Tau_lambda2)
xlabel('longueur d onde en m');
ylabel('tranmission' );
hold on;
plot(Lambda_a,Precision,'r*');
plot(Lambda_b,Precision,'r*');
hold off
grid
zoom
figure
surf (T,Longueur_d_onde,Luminance)
title('Luminance du corps noir')
ylabel('Longueur d onde en m')
xlabel('Temprature en k')
zlabel('Luminance en W/m^2/m/K/sr')
colorbar
axis ij
rotate3d
figure
surf (T,Longueur_d_onde2,Luminance_DeltaLambda)
title('Luminance du corps noir entre Lambda_a et Lmabda_b')
ylabel('Longueur d onde en m')
xlabel('Temprature en k')
zlabel('Luminance en W/m^2/m/K/sr')
colorbar
axis ij
rotate3d
figure
surf (T,Longueur_d_onde2,Luminance_DeltaLambda_Tau)
title('Luminance du corps noir entre Lambda_a et Lambda_b en fonction de \tau')
ylabel('Longueur d onde en m')
xlabel('Temprature en k')
zlabel('Luminance en W/m^2/m/K/sr')
colorbar
axis ij
rotate3d
figure
semilogy (T,Luminance_integer,T,Luminance_integer_tau,T,Luminance_Lambda,T,Luminance_integer3)
%semilogy (T,Luminance_integer3)
title('Lumiance du corps noir entre Lambda_a et Lambda_b et en fonction de \tau')
xlabel('Temprature en K')
ylabel('Luminance en W/m^2/K/sr')
legend('Luminance \lambda_a\lambda_b','Luminance \lambda_a\lambda_b \tau','Luminance centre \lambda +\delta\lambda','\int Luminance d\lambda')
grid
zoom
figure
plot(Longueur_d_onde2,Luminance_DeltaLambda(:,end),Longueur_d_onde2,Luminance_DeltaLambda_Tau(:,end))
title('Comparaison de Lumiance du corps noir entre Lambda_a et Lambda_b et en fonction de \tau pour')
xlabel('Longueur d onde en m')
ylabel('Luminance en W/m^2/m/K/sr')
legend('T=',num2str(Tmax),'\tau =',num2str(Tau))
grid
zoom
figure
semilogy (T-273,Nombre_Photon,T-273,Nombre_Photon_Lambda)
%bar(T,Nombre_Photon)
title('Nombre de photon emis par un corps noir  ')
xlabel('Temprature en C')
ylabel('Nombre de Photon en s^-1/m^2/sr ')
legend('entre \lambda_a \lambda_b \tau','centre \lambda +d\lambda')
grid
zoom

% fin de la boucle tant que
boucle = boucle + 1;
n = n + 1;

end
Resultat; % Recupration de la matrice rsultat
Nombre_Photon;

%Dtermination du temps de calcul
disp(['le temps d''execution est : ',num2str(cputime -t0),' s'])