image thumbnail

Exact Blackbodies Functions

by

 

20 Oct 2005 (Updated )

For calculating any exactely any irradiance.

planck7.m
%    Calculate the exact irradiance for  a scale of wavelength and   
%    the curves for a blackbody       
   
%%   Define the vectot temperature
%    Define the wavelength
%%   Range and increment.

%%We dertermine the differnce of luminance for  many wavelength.

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

clear all;
close all;
t0 = cputime;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%disp('Ce programme permet de calculer la luminance exacte suivant la loi de planck et  partir des courbes rduites 
%de ce mme planck sans pour autant ngliger le point fondementale de comparaison entre les longueurs d'onde fixes ainsi
% que le delta de longueur d'onde.  Pour cela on a cre deux vecteur le premier sur le delta de longueur d'onde et le second 
% sur les longueurs d'onde mme. On calcul aprs un vecteur longueurs d'onde centre sur les valeurs dsire + ou - 
% delta lambda sur deux .La taille de ce nouveau vecteur (c'est  dire sa rpartition entre le dbut et la fin de faon linaire
%doit tre imprativement de la mme taille que le vecteur 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
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; %Premire cnstante de Planck W/m2/sr
C2 = h*c/k; % Seconde constante de Planck m.K
b = (c1  / (2898^5 *(exp (c2 / 2898)-1))); %constante b pour le calcul des courbes rduites
Sigma = b * c2 / (c2/2898.8) * 1.5203 ; % Constatne de Stefan ou plutt K et si onsouhaite avoir 5.67e-8 on multiplie par pi

Tmin = 473;
Tmax = 5273;
%Tmoy = (Tmax + Tmin)/2;
DeltaT = 100;

% Dtermination des deux vecteur Temprature et Longueur d'onde et initialisation 
T = 0;  
T = Tmin: DeltaT :Tmax;
nbr = length(T);
Lambda_max = 2898.8 ./ 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 = [0.002;0.02;0.2];
B=length(D_Longueur_d_onde);
Longueur_d_onde = linspace(0.45, 3, 3)
%Longueur_d_onde = [0.36; 0.4; 0.45; 0.5; 0.55; 0.6]
A=length(Longueur_d_onde);
Longueur_d_onde1=0;

l=1;
for i=1:B
   
   for j=1:A
      
      Longueur_d_onde1 = linspace(Longueur_d_onde(j)- D_Longueur_d_onde(i)/2,Longueur_d_onde(j)+ D_Longueur_d_onde(i)/2,nbr);
      
      for k=l % cette boucle permet de remplir une matrice de longueur d'onde sans besoin de dbut ni de fin
         matrice_longueur_d_onde(:,k) = Longueur_d_onde1'; %Elle s'arrte automatiquement lorsque les boucles prcedentes
         l=l+1 ; % arrivent  leur fin
      end
      
      
   end
   
end
A=0;
B=0;
i=0;
j=0;
l=0;
matrice_longueur_d_onde;

%%par les courbes rduites
A =  c1 / b * (1 / 2898)^5;
xx = logspace(-2,2,1000); % Cration de l'abscise de la courbe rduite de planck 
x_max = length(xx);

for y = 1:x_max
   x = xx(y);
   yy(y) = A * x^(-5)/ (exp (c2/2898/x)-1); % Calcul de l'ordonne de la courbe rduite de planck
end
SumTotalY  = trapz(xx,yy); % calcul intgrale totale de la courbe rduite de planck correspond  la litrature 1.5207

SumXaYa = 0;
SumXbYb = 0;
z_a = 0;
z_b = 0;
dz = 0;
Luminance_Lambda_a_Lambda_b = 0;
d_z_b_z_a =0 ;

for i=1:k
   Lambda_a(i) = matrice_longueur_d_onde(1,i); %On prend la premire et la dernire ligne de la matrice longueur d'onde 
   Lambda_b(i) = matrice_longueur_d_onde(end,i); 
   
   for j=1:length(Lambda_max)
      
      x_a(j,i) = Lambda_a(i) / Lambda_max(j);% Pour ainsi calcul  chaques valeur la coordonnes rduites correspondantes  
      x_b(j,i) = Lambda_b(i) / Lambda_max(j);%pour toutes les valeurs du vecteur temprature
            
      xx_x_a = 0;
      xx_x_b = 0;
      yy_y_a = 0;
      yy_y_b = 0;  
            
      xx_x_a  = linspace(0.01,x_a(j,i),x_max);% Pour chaques valeur en coordonne rduites on dtermine un nouveau 
      % vecteur qui permet de calculer des intgrales plus prcisement
      for m=1:length(xx_x_a)
         x = xx_x_a(m);
         yy_y_a(m)  = A * x^(-5)/ (exp (c2/2898/x)-1);
      end
      
      x=0;
      xx_x_b  = linspace(0.01,x_b(j,i),x_max);
      
      for n=1:length(xx_x_b)
         x = xx_x_b(n);
         yy_y_b(n)  =  A * x^(-5)/ (exp (c2/2898/x)-1);
      end
      x=0;
      
      xx_x_a;
      xx_x_b;
      yy_y_a;
      yy_y_b; 
      
      SumXaYa = trapz(xx_x_a,yy_y_a);
      SumXbYb = trapz(xx_x_b,yy_y_b);      
      z_a = SumXaYa /  SumTotalY;
      z_b = SumXbYb /  SumTotalY;
      d_z_b_z_a  = z_b - z_a;     
      Luminance_totale(j) = Sigma * T(j).^4;
      Luminance_Lambda_a_Lambda_b(j,i)  = Luminance_totale(j)* d_z_b_z_a ; % Calcul final
      
   end
   
   Luminance_Lambda_a_Lambda_b;
   
end

Luminance_Lambda_a_Lambda_b

figure
plot(xx,yy,xx_x_a,yy_y_a,'g*',xx_x_b,yy_y_b,'ro')
grid
zoom
axis ([xx(1) 5 0 1]);

figure
semilogy(T,Luminance_Lambda_a_Lambda_b)
grid
zoom
legend('1')

figure
semilogy(T,Luminance_Lambda_a_Lambda_b(:,1),T,Luminance_Lambda_a_Lambda_b(:,4),T,Luminance_Lambda_a_Lambda_b(:,7))
grid
zoom
title('Calcul de la luminance pour une longueur \lambda= 0.45\mum \pm\delta\lambda\div 2')
xlabel('Temprature en K')
ylabel('Luminance en W/m^2/\mum^4/sr')
legend(num2str(D_Longueur_d_onde(1)),num2str(D_Longueur_d_onde(2)),num2str(D_Longueur_d_onde(3)))
figure
semilogy(T,Luminance_Lambda_a_Lambda_b(:,2),T,Luminance_Lambda_a_Lambda_b(:,5),T,Luminance_Lambda_a_Lambda_b(:,8))
grid
zoom
title('Calcul de la luminance pour une longueur \lambda= 1.725\mum \pm\delta\lambda\div 2')
xlabel('Temprature en K')
ylabel('Luminance en W/m^2/\mum^4/sr')
legend(num2str(D_Longueur_d_onde(1)),num2str(D_Longueur_d_onde(2)),num2str(D_Longueur_d_onde(3)))
figure
semilogy(T,Luminance_Lambda_a_Lambda_b(:,3),T,Luminance_Lambda_a_Lambda_b(:,6),T,Luminance_Lambda_a_Lambda_b(:,9))
grid
zoom
title('Calcul de la luminance pour une longueur \lambda= 3.0\mum \pm\delta\lambda\div 2')
xlabel('Temprature en K')
ylabel('Luminance en W/m^2/\mum^4/sr')
legend(num2str(D_Longueur_d_onde(1)),num2str(D_Longueur_d_onde(2)),num2str(D_Longueur_d_onde(3)))

figure
semilogy(T,Luminance_Lambda_a_Lambda_b(:,1),T,Luminance_Lambda_a_Lambda_b(:,2),T,Luminance_Lambda_a_Lambda_b(:,3))
title('Calcul de la luminance pour trois longueur d onde pour \pm\delta\lambda= 0.001\mum ')
grid
zoom
legend(num2str(Longueur_d_onde(1)),num2str(Longueur_d_onde(2)),num2str(Longueur_d_onde(3)))
xlabel('Temprature en K')
ylabel('Luminance en W/m^2/\mum^4/sr')


figure
subplot(2,2,1)
semilogy(T,Luminance_Lambda_a_Lambda_b)
subplot(2,2,2)
semilogy(T,Luminance_Lambda_a_Lambda_b(:,1),T,Luminance_Lambda_a_Lambda_b(:,4),T,Luminance_Lambda_a_Lambda_b(:,7))
subplot(2,2,4)
semilogy(T,Luminance_Lambda_a_Lambda_b(:,1),T,Luminance_Lambda_a_Lambda_b(:,2),T,Luminance_Lambda_a_Lambda_b(:,3))


disp(['le temps d''execution est : ',num2str(cputime -t0),' s'])

Contact us