• View
• Edit
• History
• Print

IllumPhys

3D modelling of the artificial sky brightness in heterogeneous environments

Martin Aubé ¹

Malcolm St-John ¹

¹Département de physique, Cégep de Sherbrooke, Québec, Canada

Last update: January 2016

Introduction to ILLUMINA

The goal of ILLUMINA is to accurately simulate the artificial multi-spectral brightness of the night sky considering the heterogeneity of the environment. We hope this model will allow for a better comprehension of this complex and non-linear phenomenon.

Entries required for the model

• Information concerning the layout of light sources along with their spectrum and angular emission functions
• Map of ground reflectance
• A numerical map of the terrain elevation
• Data concerning the atmosphere:
  • the aerosol optical depth (AOE)
  • the aerosol scattering and extinction crosssections which is dependant on the aerosol model (urban, rural, maritime, other) and relative humidity
  • the ground level atmospheric pressure

Objective of this document

The aim of this document is to give technical details concerning the equations used for ILLUMINA. The theoretical equations will be presented along with some numerical implementation.

Basic equations

We will begin with the basic theory necessary for the elaboration of ILLUMINA.

Figure 1: Rayons considérés dans le calcul

The main equation for ILLUMINA is:

{## I_{tot}=I_1 + I_{r1} + I_2 + I_{r2} ~ ~ ~ ~ (1)##}

Where {#I_1#} is the the light directly scattered toward the observer by the voxel in the line of sight (center of the observers field of view (FOV)), {#I_{r1}#} is the light reflected from the ground that is subsequently scattered toward the observer by the voxel in the line of sight , {#I_2#} is the light scattered by one voxel and subsequently scattered toward the observer by the voxel in the line of sight, et {#I_{r2}#} is the light reflected by the ground and then scattered by a voxel and finally scattered toward the observer by a voxel in the line of sight.

This equation can be reformulated to concord with the path followed by each of these intensities:

{## I_{tot}=I_{s-d-o} + I_{s-g-d-o} +I_{s-m-d-o} + I_{s-g-m-d-o} ~ ~ ~ ~ (2)##}

Where subscripts m, d, s, o, and g stands respectively for any air voxel, an air voxel in the line of sight, a voxel containing a light source, a voxel containing the observer, and a ground pixel.

Figure 2: Solid angles and voxels

To simplify the following documentation, we will ignore, for the moment, the two terms involving the second order scattering. The total energetic intensity of the light {#I_{tot}#} (Watt/stéradian) at a specific position d ({#x_d#}, {#y_d#}, {#z_d#}) in the atmosphere from an artificial light source situated at ({#x_s#}, {#y_s#}, {#z_s#}) can be explained by both the radiation arriving directly from the source ({#I_{s-d}#}) and the radiation from the source and reflected by the ground ({#I_{s-g-d}#}) at ({#x_g#}, {#y_g#}, {#z_g#}).

Figure 2.2: Purely mathematical description

Therefor:

{## I_{tot} \approx I_{s-d} + I_{s-g-d} (3)##}

Let:

{## I_{tot}=L_0 \zeta(\theta_d) T_{s-d} + \sum_g {{L_0} \zeta(180-\theta_s+\varepsilon) \Omega_{s-g} \rho(\theta_1,\theta_s,\Delta \Phi) T_{s-g-d} } (4)##}

Where {# L_{0} #} is the radiant flux of the lamp in Watt, {#\zeta(\theta)#} is the angular emission function of the lamp which allow to define the radiant intensity in Watt per steradian from the radiant flux. The angular emission function is normalized so that the integral over all solid angles of the sphere gives 1. This function is for the moment averaged on around the azimuthal axis (no azimuthal dependancy). {#T_{s-d}#} it the atmospheric transmittance on the light path separating the source and the scattering voxel while {#T_{s-g-d}=T_{g-d} T_{s-g}#} is the global atmospheric transmittance for the source-surface-voxel path. {# \Omega_{s-g} #} is the solid angle of the ground pixel as seen from the source. {#\theta_1#} is the angle between the normal line to the surface and the line joining the ground pixel and the scattering voxel. {#\Phi_1#} est is tha azimuthat angle of the same line projected on the ground surface.

{# L_0 \zeta(180-\theta_s+\varepsilon) \Omega_{s-g} #} est le flux atteignant la surface de sol (en Watt) il est donc équivalent au produit de l'éclairement sur la surface par l'aire de la surface. {#\rho(\theta_1,\theta_s,\Delta \Phi)#} est la fonction de réflectance bi-directionnelle (unités de 1/sr), c'est à dire le rapport de la luminance réfléchie (W/sr.m²) à l'angle de sortie ({#\theta_1,\Phi_1#}) sur l'éclairement (W/m²) incident sur la surface à l'angle d'entrée ({#\theta_s,\Phi_s#}). Nous définissons {#\Delta \Phi = \Phi_s - \Phi_1#} comme la différence d'angle azimuthal entre les rayons d'entrée et de sortie dans le plan de la surface réfléchissante. La luminance est définie comme le flux par unité d'angle solide par unité de surface apparente. Par conservation de l'énergie, il faut s'assurer que l'intégrale de {#\rho(\theta_1,\theta_s,\Delta \Phi) #} sur la demi-sphère soit égale à la fraction de la lumière incidente qui a été réflechie par la surface. Pour une surface parfaitement blanche on aurait par exemple que l'intégrale est égale à 1.

Pour une surface lambertienne, la luminance sortante est par définition constante. C'est à dire que le flux par unité d'angle solide par unité de surface apparente est constant peu importe la direction d'observation. En d'autres termes,

{## \frac{\rho(\theta_1,\theta_s,\Delta \Phi)}{\cos(\theta_1)} = k ~ ~ ~ ~ (5)##}

où k est une constante.

Ainsi

{## \rho(\theta_1,\theta_s,\Delta \Phi) = k \cos(\theta_1) (6)##}

La conservation de l'énergie sur la demi-sphère pour une surface lambertienne blanche donne:

{## \int_{\frac{1}{2}sphere} \rho(\theta_1,\theta_s,\Delta \Phi) d \Omega = \rho_g (7)##}

{## \int_{\theta_1=0}^{\frac{\pi}{2}} \int_{\Phi=0}^{2 \pi} k \cos(\theta_1) \sin(\theta_1) d \theta_1 d \Phi = \rho_g (8)##}

{## k \int_{\theta_1=0}^{\frac{\pi}{2}} \int_{\Phi=0}^{2 \pi} \cos(\theta_1) \sin(\theta_1) d \theta_1 d \Phi = \pi k \int_{\theta_1=0}^{\frac{\pi}{2}} \sin (2 \theta_1 ) d \theta_1 = \rho_g (9)##}

{## \pi k _{0}^{\frac{\pi}{2}}\left [ - \cos(2 \theta_1 ) \right ] = k \pi = \rho_g (10)##}

d'où {# k= \frac{\rho_g}{\pi} (11)#}

{## \rho(\theta_1,\theta_s,\Delta \Phi) = \frac {\rho_g \cos(\theta_1) }{\pi} (12)##}

où {# \rho_g #} est le facteur de réflectance hémisphérique. Cette valeur varie de 0 pour un noir absolu à 1 pour une réflexion parfaite.

Nous pouvons enfin réécrire l'équation 4 comme suit pour le cas lambertien:

{## I_{tot}=L_0 \zeta(\theta_d) T_{s-d} + \frac {1}{\pi} \sum_g {{L_0} \zeta(180-\theta_s+\varepsilon) \Omega_{s-g} \rho_g \cos(\theta_r + \varepsilon) T_{s-g-d} } (13)##}

car {# \theta_1 = \theta_r + \varepsilon#}.

Flux traversant le voxel

Comme chaque voxel couvre un certain angle solide qui varie selon l'origine de l'émission de lumière, le concept de flux énergétique est plus représentatif du rayonnement traversant le voxel. Le flux énergétique est la puissance radiative traversant le voxel (W).

{## \Phi_{tot}=L_0 \zeta(\theta_d) T_{s-d} \Omega_{s-d} + \frac {1}{\pi} \sum_c {{L_0} \zeta(180-\theta_s+\varepsilon) \Omega_{s-g} \rho_c \cos(\theta_r + \varepsilon) T_{s-g-d} \Omega_{g-d} } (14)##}

{#\Omega_{s-d}#} est l'angle solide couvert par le voxel atmosphérique tel que vu de la position de la source alors que {#\Omega_{g-d}#} est l'angle solide couvert par le même voxel telle que vu de la position de la surface réfléchissante. Pour calculer ces angles solides, nous utilisons la méthodologie présentée plus bas en considérant la plus grande valeur d'angle solide parmi 3 sections du voxel selon les plans cartésiens du système de coordonnée et passant par le centre du voxel.

Figure 3: Description of voxel

Radiance au capteur

Comme nous sommes intéresses à modéliser l'intensité énergétique diffusée (W/sr) dans une direction de visée donnée nous devons multiplier chacun des termes par la probabilité de diffusion ({#P#}) par unité d'angle solide dans cette direction. Cette probabilité est donnée par le produit de la probabilité de diffusion dans le voxel avec la fonction de phase de diffusion {#F(\pi-\theta_d+\theta_v)#}. La fonction de phase doit être exprimée en unité de 1/sr, son intégration sur la sphère doit donner 1 (c.-à-d. {# \int_{2 \pi} F(\pi-\theta_d+\theta_v) d\Omega = 1 #}).

{## P = 1 - T_{voxel} (16)##}

{# T_{voxel} #} est la transmittance du voxel selon la direction de la lumière incidente. La méthode utilisée pour déterminer cette transmittance sera détaillée plus loin.

L'intensité énergétique reçue par l'observateur est simplement l'intensité énergétique diffusée réduite par la transmittance entre le voxel et l'observateur.

{## I_{obs} = I_{diffusée} T_{d-o} \gamma_{capteur} ~ ~ ~ ~ (17) ##}

On a multiplié l'intensité attribuable au voxel par la fraction du voxel entrant dans le champ de vision du capteur {#\gamma_{capteur}#}.

{## \gamma_{capteur} = \frac {\Omega_{FOV}}{\Omega_{o-d}} ~ ~ ~ ~ (18)##}

Pour obtenir le flux mesuré au capteur en provenance du voxel, il faut multiplier {#I_{obs}#} par l'angle solide couvert par le capteur {#\Omega_{capteur}#} tel que perçu du voxel diffusant.

{## \Phi_{obs} = I_{obs} \Omega_{capteur} ~ ~ ~ ~ (19)##}

L'éclairement énergétique {#E_{obs}#} (W/m²) peut être calculé en divisant la flux au capteur {#\Phi_{obs}#} par la surface du capteur {#A_{capteur}#}

{## E_{obs} = \frac {\Phi_{obs}}{A_{capteur}} ~ ~ ~ ~ (20)##}

Finalement si on s'intéresse à la radiance mesurée, il s'agit de diviser {# E_{obs}#} par l'angle solide du champ de vision du capteur {#\Omega_{FOV}#}.

{## R_{obs} = \frac {E_{obs}}{\Omega_{FOV} } ~ ~ ~ ~ (21)##}

Finalement on peut exprimer la radiance en fonction de l'intensité diffusée de l'équation 15.

Pour calculer la radiance totale au capteur, il faut faire la somme des {#R_{obs}#} de tous les voxels croisés par la ligne de visée du capteur.

Intégration de la diffusion de second ordre

Pour prendre en considération la diffusion de 2e ordre on ajoute simplement une diffusion vers un voxel intermédiaire avant d'atteindre le voxel de la ligne de visée. Le nombre de calculs requis croît très rapidement. Nous définissons un rayon d'action maximal de la 2e diffusion autour du parcours optique entre la source et le voxel de la ligne de visée. Lorsque le modèle est en configuration de calcul détaillé, tous les voxels situés dans un rayon inférieur du rayon d'action maximal de tout point de la ligne source-voxel de la ligne de visée seront calculés. Comme ce calcul est très coûteux en temps de calcul nous avons déterminé un saut de diffusion qui permet de ne réaliser qu'un calcul sur N, où N est idéalement un nombre premier.

Calcul de l'angle solide formé par la surface de sol telle que vue de la source {# \Omega_{s-g} #}

L'angle solide formé par la surface au sol à partir de la source est dépendant de la direction source-surface ainsi que de l'inclinaison de cette surface. C'est la surface apparente {# S' #} vue de la source qui est déterminante.

{## \Omega_{s-g} = \frac {S'}{{r_s}^2} ~ ~ ~ ~ (23)##}

{# r_s #} est la distance source-surface.

La surface apparente peut être calculée en prenant le produit scalaire du vecteur {# \vec S #} avec le vecteur {# \vec r_s #}. Toutefois nous avons estimé que cette façon de calculer est gourmande en temps de calcul, C'est pourquoi nous proposons une autre approche qui consiste à utiliser des identités trigonométriques sur des triangles sphériques (triangles dessinés sur la surface d'une sphère). La première étape est de déterminer l'intersection des vecteurs reliant la source et chacun des coins de la surface réfléchissante avec une sphère unitaire (voir figure ci-dessous). Cette méthode sera utilisée pour calculer l'ensemble des angles solides dans le modèle.

Figure 4: Description of solid angle

Soient les vecteurs {# \vec r_1 #}, {# \vec r_2 #}, {# \vec r_3 #}, {# \vec r_4 #} menant à chacun des coins de la surface données par :

{## \vec r_1 = \left ( x_g-\frac{\Delta x}{2} -x_s \right ) \vec{i} + \left ( y_g+\frac{\Delta x}{2} -y_s \right ) \vec{j} + \left ( z_g - \tan(\varepsilon_i) \frac{\Delta x}{2} + \tan(\varepsilon_j) \frac{\Delta x}{2} -z_s \right ) \vec{k} ~ ~ ~ ~ (24)##}

{## \vec r_2 = \left ( x_g+\frac{\Delta x}{2} -x_s \right ) \vec{i} + \left ( y_g+\frac{\Delta x}{2} -y_s \right ) \vec{j} + \left ( z_g + \tan(\varepsilon_i) \frac{\Delta x}{2} + \tan(\varepsilon_j) \frac{\Delta x}{2} -z_s \right ) \vec{k} ~ ~ ~ ~ (25)##}

{## \vec r_3 = \left ( x_g-\frac{\Delta x}{2} -x_s \right ) \vec{i} + \left ( y_g-\frac{\Delta x}{2} -y_s \right ) \vec{j} + \left ( z_g - \tan(\varepsilon_i) \frac{\Delta x}{2} - \tan(\varepsilon_j) \frac{\Delta x}{2} -z_s \right ) \vec{k} ~ ~ ~ ~ (26)##}

{## \vec r_4 = \left ( x_g+\frac{\Delta x}{2} -x_s \right ) \vec{i} + \left ( y_g-\frac{\Delta x}{2} -y_s \right ) \vec{j} + \left ( z_g + \tan(\varepsilon_i) \frac{\Delta x}{2} - \tan(\varepsilon_j) \frac{\Delta x}{2} -z_s \right ) \vec{k} ~ ~ ~ ~ (27)##}

Nous considérons les triangles sphériques formés par l'intersection avec la sphère des triplets de vecteurs {# \vec r_1 #}, {# \vec r_2 #}, {# \vec r_3 #} puis {# \vec r_2 #}, {# \vec r_3 #}, {# \vec r_4 #}. Les cotés du premier triangle sont donnés par

{## a = \arccos{\left(\frac{\vec r_2\cdot\vec r_3}{r_2 r_3 } \right)} ~ ~ ~ ~ (28)##}

{## b = \arccos{\left(\frac{\vec r_1\cdot\vec r_3}{r_1 r_3 } \right)} ~ ~ ~ ~ (29)##}

{## c = \arccos{\left(\frac{\vec r_1\cdot\vec r_2}{r_1 r_2 } \right)} ~ ~ ~ ~ (30)##}

Si nous posons,

{## s=\frac{a+b+c}{2} ~ ~ ~ ~ (31)##}

le théorème de Huilier nous permet de calculer la surface de ce triangle:

{## A_{123}= 4 \tan^{-1} \left(\sqrt{\tan\left(\frac{s}{2}\right) \tan\left(\frac{s-a}{2}\right) \tan\left(\frac{s-b}{2}\right) \tan\left(\frac{s-c}{2}\right)} \right) ~ ~ ~ ~ (32)##}

Pour le second triangle, nous obtenons un résultat similaire:

{## a = \arccos{\left(\frac{\vec r_2\cdot\vec r_3}{r_2 r_3 } \right)} ~ ~ ~ ~ (33)##}

{## d = \arccos{\left(\frac{\vec r_3\cdot\vec r_4}{r_3 r_4 } \right) } ~ ~ ~ ~ (34)##}

{## e = \arccos{\left(\frac{\vec r_2\cdot\vec r_4}{r_2 r_4 } \right)} ~ ~ ~ ~ (35)##}

{## s'=\frac{a+d+e}{2} ~ ~ ~ ~ (36) ##}

{## A_{234}= 4 \tan^{-1} \left(\sqrt{\tan\left(\frac{s'}{2}\right) \tan\left(\frac{s'-a}{2}\right) \tan\left(\frac{s'-d}{2}\right) \tan\left(\frac{s'-e}{2}\right)} \right) ~ ~ ~ ~ (37)##}

La somme de la surface des deux triangles donne la surface de la portion de sphère bornée par les 4 intersections. Comme la sphère possède un rayon unitaire, cette surface est égale à l'angle solide apparent de la surface S.

{## \Omega_{s-g} = A_{123} + A_{234} ~ ~ ~ ~ (38)##}

Calcul de l'aire de la surface réfléchissante inclinée {#\vec S #}

Soit la surface {#\vec S#} définie par le produit scalaire

{## \vec S = \vec m_i \times \vec m_j ~ ~ ~ ~ (39)##}

avec

{## \vec m_i = \Delta x (\vec i + \tan(\varepsilon_i) \vec k) ~ ~ ~ ~ (40)##}

et

{## \vec m_j = \Delta x (\vec j + tan(\varepsilon_j) \vec k) ~ ~ ~ ~ (41)##}

où {#\Delta x#} est la largeur horizontale d'un voxel du modèle. {#\varepsilon_i#} est l'inclinaison de la surface par rapport à l'axe des x alors que {#\varepsilon_i#} est l'inclinaison par rapport à y.

On obtiens donc

{## \vec S = {\Delta x}^2 (-\tan(\varepsilon_i) \vec i - \tan(\varepsilon_j) \vec j + \vec k) ~ ~ ~ ~ (42)##}

puis

{## S = {\Delta x}^2 \sqrt{ \tan^2(\varepsilon_i) + \tan^2(\varepsilon_j) + 1} ~ ~ ~ ~ (43) ##}

Figure 5: Calculation of the air of an inclined surface

Calcul de l'inclinaison de la surface par rapport à l'axe z ({#\varepsilon#})

Par définition,

{## \vec k \cdot \vec S = \left | \vec S \right | \cos(\varepsilon) ~ ~ ~ ~ (44)##}

{## \cos(\varepsilon) = \frac{1}{\sqrt{\tan^2(\varepsilon_x) + \tan^2(\varepsilon_y) +1 }} ~ ~ ~ ~ (45)##}

{## \varepsilon = \arccos{\left(\frac{1}{\sqrt{\tan^2(\varepsilon_x) + \tan^2(\varepsilon_y) +1 }} \right)} ~ ~ ~ ~ (46)##}

Calcul de l'aire apparente de la surface telle que perçue du voxel ({#\vec S'#})

Le calcul de cette aire est important car le rapport {# \frac{S'}{S} = \cos(\theta_r + \varepsilon) #} intervient dans l'équation 13.

{## S' = \frac{1}{r_d} \vec S \cdot \vec r_d ~ ~ ~ ~ (47)##}

{## \vec r_d = (x_d-x_g) \vec i + (y_d-y_g) \vec j + (z_d + H_s - z_g) \vec k ~ ~ ~ ~ (48)##}

d'où

{## r_d = \sqrt{~(x_d-x_g)^2+(y_d-y_g)^2+(z_d+H_s-z_g)^2} ~ ~ ~ ~ (49) ##}

En substituant les équations 43, 48 et 49 dans l'équation 47, nous obtenons

{## S' = {\Delta x}^2 \left( \frac {-(x_d-x_g) \tan(\varepsilon_i) -(y_d-y_g) \tan(\varepsilon_j) + (z_d + H_n - z_g)}{\sqrt{(x_d-x_g)^2+(y_d-y_g)^2+(z_d + H_s - z_g)^2}} \right) ~ ~ ~ ~ (50)##}

En combinant l'équation 43 avec la formulation de {#S'#} ci-dessus, nous obtenons:

{## \frac{S'}{S} = \cos(\theta_r + \varepsilon)= \frac{-(x_d-x_g) \tan(\varepsilon_i) -(y_d-y_g) \tan(\varepsilon_j) + (z_d + H_s - z_g)}{\sqrt{(x_d-x_g)^2+(y_d-y_g)^2+(z_d + H_s - z_g)^2}\sqrt{\tan^2(\varepsilon_i) + \tan^2(\varepsilon_j) + 1} } ~ ~ ~ ~ (51) ##}

Normalisation de la fonction d'émission

Comme la fonction {#\zeta(\theta)#} sert à passer d'une luminosité en Watt à une intensité énergétique en Watt/sr il est essentiel que cette fonction soit normalisée sur la sphère de sorte que son intégrale sur toutes les directions soit égale à 1 car:

{## \int_0^\pi L_0 \zeta(\theta) d\Omega = L_0 ~ ~ ~ ~ (52)##}

{## \int_0^\pi \zeta(\theta) d\Omega = 1 ~ ~ ~ ~ (53)##}

Avec, {## d\Omega = \frac{2 \pi r \sin(\theta) r d\theta}{r^2} = 2 \pi \sin(\theta) d\theta ~ ~ ~ ~ (54)##}

Figure 6: Solid angle for a band around a sphere

Calcul de la transmittance

La transmittance totale le long d'un parcours ({#T#}) est obtenue par le produit des transmittances liées à la diffusion par les molécules gazeuses {#T_R#} (diffusion de Rayleigh), l'absorption par les molécules gazeuses {#T_A#} et l'atténuation par les aérosols (diffusion et absorption) {#T_a#}.

{## T=T_R T_A T_a ~ ~ ~ ~ (55)##}

Comme le modèle est discret (voxels en 3D), il s'avère plus pratique d'exprimer la transmittance totale le long du parcours comme le produit des transmittances de chaque couches individuelles.

{## T=T_1 T_2 T_1 ... T_N ~ ~ ~ ~ (56)##}

Diffusion de Rayleigh

Nous supposons que le profil vertical de densité {#n#} et de pression {#p#} des molécules stables épousent une décroissance exponentielle avec l'altitude selon une échelle de hauteur {#H_R=8000m#}.

{## \frac{n}{n_0} = \frac{p}{p_0} = e^{-z/H_R)} ~ ~ ~ ~ (57)##}

Soit la formulation de {#T_R#} donnée par Kneizys et al. (1980),

{## T_{R\infty} = \exp\left( \frac {-M'}{\lambda^4 \left( 115.6406- \frac{1.335}{\lambda^2} \right)} \right) ~ ~ ~ ~ (58)##}

i.e. {## T_{R\infty} = \exp\left( -M' \tau_R \right) ~ ~ ~ ~ (59) ##}

où {## \tau_R = \frac {1}{\lambda^4 \left( 115.6406- \frac{1.335}{\lambda^2} \right)} ~ ~ ~ ~ (60)##}

{#M'#} est la masse d'air corrigée pour la variation de pression au sol {#p'_0#} par rapport à la pression normale {#p_0#}

{## M' = M \frac{p'_0}{p_0} ~ ~ ~ ~ (61)##}

La masse d'air {#M#} est donnée par Kasten (1966)

{## M=\frac{1}{\cos\theta + 0.15 (93.885-\theta)^{-1.253}} ~ ~ ~ ~ (62)##}

La masse d'air est proportionnelle à l'intégrale de la densité du gaz le long du parcours optique et qu'elle est normalisée de telle sorte que la masse d'air totale en visée verticale vaut 1.

Si nous appliquons cette idée à un voxel du niveau vertical k, nous obtenons

{## \Delta M_k = \frac{\int_{z_k}^{z_{k+1}} \exp \left(-z/H_R \right) \frac{1}{\cos(\theta)}dz}{\int_{0}^{\infty} \exp \left(-z/H_R \right) dz} (63)##}

ou

{## \Delta M_k = \frac{1}{\cos(\theta)} \frac{\int_{z_k}^{z_{k+1}} \exp \left(-z/H_R \right)dz}{\int_{0}^{\infty} \exp \left(-z/H_R \right) dz} ~ ~ ~ ~ (64)##}

Cette équation peut être généralisée à une plus grande épaisseur d'atmosphère en remplaçant {#z_{k+1}#} par {#z_{k'}#}.

Pour une visée à travers l'atmosphère au complet nous obtenons :

{## M = \frac{1}{\cos(\theta)} \frac{\int_{0}^{\infty} \exp \left(-z/H_R \right)dz}{\int_{0}^{\infty} \exp \left(-z/H_R \right)dz} = \frac{1}{\cos(\theta)} ~ ~ ~ ~ (65)##}

Il s'agit de l'équation approchée de la masse d'air (atmosphère plan parallèle). Pour corriger pour la courbure de la Terre, il faut simplement remplacer le facteur {#\frac{1}{\cos(\theta)}#} par l'équation 62. Toutefois comme nous appliquons cette méthodologie à l'échelle régionale, nous pouvons faire l'approximation de l'atmosphère plan parallèle.

Un problème survient lorsque {#\theta=90^o#} car {#1/\cos\theta#} tend vers l'infini. Nous corrigeons ce problème en traitant ce cas particulier en remplaçant l'intégration verticale par une intégration horizontale.

{## \Delta M_k = \frac{\int_{x_k}^{x_{k+1}} \exp \left(-z_k/H_R \right)dx}{\int_{0}^{\infty} \exp \left(-z/H_R \right)dz} (66)##}

{## \Delta M_k = \frac{|x_k-x_{k+1}| \exp \left(-z_k/H_R \right)}{\int_{0}^{\infty} \exp \left(-z/H_R \right)dz} (67)##}

Étant donné que

{## \int_{0}^{\infty} \exp \left(-z/H_R \right) dz = H_R (68)##}

et que

{## \int_{z_k}^{z_{k+1}} \exp \left(z/H_R \right) dz = H_R \left( \exp \left(-z_k/H_R \right)- \exp \left(-z_{k+1}/H_R \right)\right) (69)##}

Nous pouvons donc exprimer la transmittance du parcours optique k comme:

{## T_{Rk} = \exp(-\tau_R \Delta M'_k ) (72)##}

où l'élément de masse d'air corrigé pour les variations de pression atmosphérique au sol est donné par:

{## \Delta M'_k = \Delta M_k \left(\frac{p'_0}{p_0} \right) (73)##}

Diffusion par les aérosols

Pour le cas des aérosols nous disposons de l'épaisseur optique des aérosols pour l'atmosphère au complet {#\tau_a#}. Nous pouvons appliquer l'équation 72 en remplaçant {#\tau_R#} par {#\tau_a#} pour calculer la transmittance liée aux aérosols. De cette transmittance nous pouvons déterminer la probabilité de diffusion par les aérosols (équation 16). La probabilité de diffusion est alors la somme de la probabilité de diffusion par les aérosols et la probabilité de diffusion par les molécules.

Fonctions de phase de diffusion

Fonction de phase des molécules

Fonction de phase des aérosols

Absorption par les gaz

Nous choisissons les longueurs d'ondes de modélisation de façon à éviter les régions du spectre où se produit l'absorption moléculaire afin de ne pas avoir à modéliser ce processus relativement complexe. Le choix des longueurs d'ondes optimales est illustré par les zones ombragées sur la figure ci-dessous.

Cloud contribution to the total radiance

The goal of this work is to create a sub routine for the light pollution model ILLUMINA. This ad on will allow the users to model the effect of clouds on light pollution both above and below clouds.

Sky brightness below clouds

To incorporate clouds to the already existing model, it would suffice to enter the height of the clouds as the maximum height of the atmosphere and add a reflectance function for the cloud base according to the cloud type. The reflectance of different types of clouds has already been studied by Shapiro (1982).

Table 1: Estimates of reflectivity for overcast cloud layers

Cloud Type	{#\cos(z)#}	0.05	0.15	0.25	0.35	0.45	0.55	0.65	0.75	0.85	0.95
Thin Ci/Cs R1		0.240	0.235	0.211	0.170	0.143	0.134	0.120	0.106	0.103	0.100
Thick Ci/Cs R1		0.568	0.535	0.460	0.401	0.342	0.297	0.273	0.249	0.240	0.240
As/Ac R2		0.650	0.645	0.621	0.604	0.587	0.576	0.569	0.562	0.560	0.560
Cu/Cb R3		0.661	0.650	0.636	0.611	0.592	0.575	0.561	0.549	0.535	0.520
Sc/St R3		0.703	0.693	0.676	0.660	0.648	0.639	0.629	0.620	0.613	0.609

Where z is the incident zenith angle of the light to the clouds. R is the reflectance and the number following it is the layer of the cloud as given by table 2. The clouds are listed as follows: Cirrus (Ci), Cirrostratus (Cs), Altostratus (As), Altocumulus (Ac), Cumulus (Cu), Cumulonimbus (Cb), Stratus (St), Stratocumulus (Sc).

Table 2: Height of cloud base

Level	Cloud type	Cloud base (km)
High	Ci/Cs R1	9.3
Middle	As/AC R2	4.0
Low	Cu/Cb R3	1.2
Low	Sc/St R3	1.1

In table 1, {#\cos(z)#} increases by bounds of 0.10. The software Fityk was used to estimate a half Gaussian curve that would allow us to have a continuous function to be used for any incident angle z (between 0 and {#\frac{\pi}{2}#}).

Table 3: Parameters of the fitted fonction to the table 1 reflectance data

Cloud Type	Constant ({#A#})	Height ({#a_0#})	Center ({#a_1#})	HWHM ({#a_2#})
Thin Ci/Cs R1	0.1	0.146	0	0.357
Thick Ci/Cs R1	0.24	0.3327	0	0.344
As/Ac R2	0.56	0.0932	0	0.34
Cu/Cb R3	0.515	0.144	0	0.49
Sc/St R3	0.61	0.0907	0	0.4087

The form of the function would therefor be:

{##\rho_{cloud}=A+ a_0 \exp\left[ - \left(\frac{\cos(z) - a_1}{a_2} \right)^2 \right] ##}

Where {#A#} is the constant, {#a_0#} is the height, {#a_1#} is the center and is {#a_2#} the half width half maximum (HWHM).

Seeing as the center is zero for each of the curves, the equation becomes :

{##\rho_{cloud}=A+ a_0 \exp\left[ - \left(\frac{\cos(z)}{a_2} \right)^2 \right] ##}

Sky brightness above clouds

TBD

Références

Kasten, F., A New Table and Approximate Formula for Relative Optical Air Mass, Arch. Meteorol. Geophys. Biochlimatol., Ser. B14, 1966, pp. 206-223.

Kneizys, F. X., E. P. Shettle, W. O. Gallery, J. H. Chetwynd, Jr., L. W. Abrea, J. E. A. Selby, R. W. Fenn, and R. W. McClatchey, Atmospheric Transmittance/Radiance: Computer Code LOWTRAN5, Tech. Rep. AFGL-TR-800067, Bedford, MA: U.5. Air Force Geophysics Laboratory, 1980.