Diffusion de photons dans un nuage plan / diagrammes de phase

Définition: Photon : Le photon est la particule élémentaire de pure énergie (sans masse). La lumière, les ondes électromagnétiques sont toutes constituées de photons.

Diffusion : Pour un photon donné, une diffusion est une collision avec un objet (poussière) qui le dévie de sa trajectoire rectiligne initiale.

Diffusion isotrope : C'est le cas où la direction du vecteur d'onde du photon est répartie uniformément dans toutes les directions de l'espace. On néglige ici le phénomène d'absorption.

Diffusion Rayleigh : Mode de diffusion où la longueur d'onde est très supérieure à la taille des particules diffusantes (diffusion élastique sans variation d'énergie).

Fluctuations statistiques : L'utilisation d'une méthode de Monte Carlo implique l'utilisation de variables aléatoires, en l'occurence les angles de diffusion et les distances parcourues des photons. La nature aléatoire de ces variables implique un bruit statistique.

Libre parcours moyen : Distance moyenne que parcourt une particule entre deux collisions.

Fonction de répartition : La fonction de répartition d'une variable aléatoire réelle caractérise la loi de probabilité de cette variable aléatoire réelle.

Générateur de nombre aléatoire : Suite pseudo-aléatoire générée par l'ordinateur qui permet de tirer des nombres aléatoires de manière uniforme.

Liste des paramètres d'entrée de l'applet

Liste des paramètres de sortie de l'applet

Mode d'emploi de l'applet

Explications

Principe

Nous allons décrire dans cette section l'algorithme qui nous a permis de réaliser cette statistique sur les trajectoires de sortie des photons.
On considère un faisceau de photons dirigé vers un nuage diffusif représenté dans notre cas par une couche à deux dimensions de largeur fixée par l'utilisateur. On exprime la largeur de la couche diffusive en unité de libre parcours moyen (λ), au lieu d'utiliser la densité de particules et la section efficace. Cette approche permet de simplifier le problème.

Schéma explicatif

La trajectoire du photon est déterminée via deux paramètres : une distance et un couple d'angles qui sont tous tirés aléatoirements grâce à un générateur de nombre pseudo-aléatoires. La position est ainsi connue au fur et à mesure des diffusions successives. Lorsque le photon sort du nuage (sa position par rapport à l'axe z est au-delà de la largeur du nuage), on passe au photon suivant.

La distance parcourue par le photon est tirée aléatoirement suivant une loi de Poisson (loi exponentielle). Pour cela on s'intéresse à la fonction de répartition telle que : f(x)->F(u).
Si l'on considère F comme étant la fonction de répartition d'une variable aléatoire continue X, alors la variable aléatoire F(X) est uniforme entre 0 et 1. On sait que si deux fonctions de répartition sont identiques, alors les variables aléatoires associées sont égales. On peut donc écrire $X=F^{-1}(U)$ . On peut utiliser un changement de variable pour générer des nombres aléatoires X sous une loi F. On part de la loi exponentielle avec F(X)=1-exp{(-x)}

. On peut en déduire :

$U=1-e^{-X}$

$e^{-X}=1-U$

$\rightarrow$ U étant tiré au hasard par le générateur, X est donc un nombre aléatoire entre 0 et 1.

Si la distance tirée est supérieure à la largeur du nuage, le photon n'est pas diffusé (il traverse directement le nuage sans interaction) et on passe alors au photon suivant.
Si la distance tirée est inférieure à la largeur du nuage (c'est à dire dans le nuage), le photon est diffusé. On tire alors un angle déterminant la direction post-diffusion du photon.

Remarque : Afin d'optimiser le temps de calcul, nous avons utilisé les probabilités conditionnelles afin qu'un photon lancé soit un photon diffusé. Cela permet de gagner un temps considérable au niveau du calcul. On exprime cette probabilité ainsi :
$P( X | D)=U=\frac{P(X \cap D)}{P(D)}$
On en déduit une loi de tirage pour la première distance de pénétration dans le nuage :
X=-log(1-U(1-exp(-D)))

On s'assure ainsi que tous les photons soient diffusés lorsqu'ils entrent dans le nuage diffusif.

La distance étant tirée il faut à présent générer deux angles (θ, phi2

) qui détermineront la direction du photon après sa diffusion.

Afin de suivre cette trajectoire, nous considérons deux référentiels distincts : le référentiel du point d'entrée dans le nuage (voir schéma 1) ainsi qu'un référentiel centré sur la particule diffusive.

Le référentiel d'entrée dans le nuage (R): il possède une base orthonormée directe (e₁,e₂,e₃), avec e₃ (axe z sur le schéma 1) perpendiculaire au plan formé par le nuage. C'est par rapport à ce repère que l'on positionne le photon dans le nuage.

Le référentiel de diffusion (R'): ce référentiel possède lui aussi une base orthonormée directe (e₁',e₂',e₃') avec e₃' dirigé suivant le vecteur d'onde du photon avant diffusion.

Si on considère les angles (θ, phi2

) comme étant les coordonnées sphériques du vecteur d'onde unitaire dans R, on peut déterminer entièrement le référentiel R' suivant deux rotations successives.
La génération aléatoire de ces deux angles fait intervenir le type de la diffusion. Dans cet applet, nous traitons deux cas seulement (bien qu'il en existe d'autres): la diffusion isotrope et la diffusion Rayleigh. La diffusion isotrope, comme son nom l'indique, produit des angles sans directions privilégiées. La diffusion Rayleigh, quant à elle, a tendance à focaliser les angles de diffusion sur les directions "avant" et "arrière" du vecteur d'onde du photon avant diffusion. Les deux prochaines sections sont dédiées à la description des algorithmes permettant de générer ces angles dans chacun des types de diffusion abordés ici.

Diffusion Isotrope

Dans le cas de la diffusion isotrope, du fait de la symétrie, la seule variable angulaire que nous avons besoin de générer est $\theta$ .
Pour θ, la fonction de répartition est:

$F(\theta) = \int_{0}^{\theta} \frac{1}{2} \sin t dt = \frac{1}{2}(1 - \cos \theta)$

Si U est une variable aléatoire uniforme entre [0,1]

, alors:
$U = \frac{1}{2}(1-\cos\theta)$

$\cos \theta =1-2U$

On en déduit que $cos(\theta)$ est donc uniforme entre -1 et 1. Par conséquent, nous utiliserons cette variable plutôt que $\theta$ . Il y a par contre un désavantage à cause de la présence des pôles en $\theta=0$ et $\theta=\pi$ . Pour s'en affranchir, nous plongeons la sphère dans l'espace euclidien $\Re^3$ , et l'on rammène tous les points tirés aléatoirement sur la sphère unité.

Comme nous l'avons vu, on tire de nouveau une distance, et on vérifie la position relative au nuage du photon. Si le photon est toujours dans le nuage, on retire un angle. Sinon on passe au photon suivant jusqu'à ce que tous les photons soient sortis du nuage.

Comme nous l'avons expliqué précédemment, nous choisissons que le photon entre dans le nuage de manière perpendiculaire, c'est à dire θ = 0 et donc cos θ=1. De cette façon, l'angle du photon est confondu avec l'axe z (cf : schéma 1).

Il est affiché à l'écran le taux de diffusion c'est à dire le rapport entre photons diffusés et photons lancés, le nombre moyen de diffusion que subit un photon, l'avancement du calcul et sa durée.

Le but de cette applet est de tracer le diagramme de diffusion. Pour cela nous devons enregistrer les coordonnées $(\theta,\phi)$ du vecteur d'onde sortant du nuage.

Le diagramme de phase est la représentation graphique de l'ensemble des angles $(\theta,\phi)$ de sortie du nuage des photons. Pour tracer ce diagramme, on échantillonne $\theta$ sur un certain nombre d'intervalles. De nos angles de sortie, on en déduit quel intervalle incrémenter.

On trouve l'indice de l'intervalle de l'angle de sortie avec la formule suivante :
$indice=\frac{cos\theta + 1}{\delta}$

On incrémente le nombre de photons pour cet intervalle. Afin de pouvoir afficher, on normalise l'histogramme au nombre de diffusion moyen par intervalle.

Dans un cas idéal où le nuage est très fin, le diagramme de phase tend vers un cercle. Vous avez la possibilité de tracer la figure théorique (le cercle) et le diagramme de phase sur la même figure (voir capture 1).

Capture d'écran 1

Diffusion Rayleigh

Dans le cas de la diffusion Rayleigh, la fonction de phase est : $f(\alpha)=\frac{3sin(\alpha)(1 +cos(\alpha)^2)}{8}$

On intègre la fonction de phase et on obtient la fonction de répartition : $F(\alpha)=\frac{4-cos(\alpha)(3 +cos(\alpha)^2)}{8}$

Il se trouve que cette fonction est relativement proche de la droite $y=\pi \alpha$ . On peut donc inverser la fonction de répartition (calculée à partir de la fonction de phase) en utilisant une méthode de Newton. On génère un nombre aléatoire u entre 0 et 1. On en déduit un angle $\alpha=\pi u$ compris entre $[0,\pi]$ . L'élément de correction $d \alpha$ est de la forme:

$d \alpha = \frac{F - u}{f}$

Puis on fait converger la valeur de l'angle $\alpha$ afin que la distribution statistique globale de ces angles suive bien une loi de diffusion de Rayleigh:

$\alpha = \alpha - \epsilon \times d \alpha$

Avec $\epsilon$ un critère de convergence fixé à 0.8 dans notre cas.
Pour obtenir les angles $(\theta,\phi)$ , on fait subir une rotation au couple $(\alpha,\beta)$ (avec $\beta$ un angle généré aléatoirement entre $[0,2\pi]$ ) conformément à ce qui a été décrit plus haut. Cette rotation nous permet de passer du référentiel du photon au référentiel d'entrée.
Enfin on retire une distance suivant la loi de Poisson, si le photon est toujours dans le nuage on réalise la même méthode expliquée ci-dessus, sinon un autre photon est envoyé. Le processus de création du diagramme de phase est strictement identique à celui décrit dans la section diffusion isotrope.
Tout comme dans le cas de la diffusion isotrope, il est possible de superposer la fonction de phase au diagramme de phase obtenu par le calcul (voir capture 2).

Capture d'écran 2

Parallélisation

Pour optimiser le temps de calcul, une parallélisation du code a été effectuée. Elle permet donc de partager le calcul entre les différents processeurs de la machine en créant des sous-processus (threads). Un sous-processus est un fil d'instructions et d'exécution (voir schéma 2).

schéma 2

Remarque : Un processeur est une puce qui contient entre autres une unité de calcul (CPU). Il peut y avoir un ou plusieurs CPU sur le même processeur. Ces différents CPU sont appelés les coeurs du processeur.

A la compilation, sans instruction particulière dans le code source, les tâches sont effectuées par un seul coeur. Si le code est parrallélisé, on peut affilier un sous-processus à chacun des coeurs. Néanmoins il est nécessaire d'effectuer des synchronisations à intervalles de temps réguliers afin de savoir où en est le calcul. Une synchronisation est un moyen de regrouper les différents résultats de tous les sous-processus, et de les envoyer vers l'entité principale de calcul. Cela correspond simplement à un temps durant lequel un sous-processus qui a fini son calcul attend que les autres aient terminé leurs calculs respectifs.

Exemple

L'applet a été exécutée sur une machine possédant un processeur 4 coeurs. Nous avons comparé dans le tableau suivant le temps de calcul pour une série de paramètres donnés:

nombre de photons diffusés : 2E7
largeur du milieu diffusif : 0.01
nombre d'intervalles pour le diagramme de phase : 1000
temps entre chaque affichage : 1 seconde

nombre de sous-processus (coeurs utilisés)	1	2	3	4
temps de calcul (secondes)	61.32	32.51	25.07	17.99

Le cas zéro sous-processus

Il est possible dans l'applet de définir zéro sous-processus. Cela signifie simplement qu'aucun sous-processus n'est créé pour réaliser le calcul. Cela permet nottament de s'affranchir du temps d'attente entre processus, et donc de gagner un très léger temps de calcul par rapport au cas à un sous-processus.

Diffusion de photons dans un nuage plan / diagrammes de phase

Table des matières

Introduction

Rappels de cours

Définition