Diffraction de la lumière

2.a. Lumière monochromatique : théorie scalaire

Une onde électromagnétique est décrite par un champ électromagnétique $(\overrightarrow{E},\overrightarrow{B})$. Mis à part l'écran diffractant, on suppose que l'onde se propage dans un milieu linéaire isotrope et homogène, d'indice de réfraction n. L'étude est menée pour une onde monochromatique de fréquence $\nu$. Chacune des 6 composantes du champ électromagnétique est alors une fonction de la position (point P) et du temps de la forme suivante :

$$u(P,t)=A\left(P\right)cos(\phi \left(P\right)-2\pi \nu t)\text{(1)}$$

Chacune de ces composantes obéit à l'équation des ondes suivante :

$${\nabla }^{2}u-\frac{n^2}{c^2}\frac{{\partial }^{2}u}{\partial t^2}=0\text{(2)}$$

Pour un milieu non dispersif, comme l'air, l'indice n est une constante. Pour un milieu dispersif, l'indice de réfraction dépend de la fréquence.

Les six composantes du champ électromagnétique se propagent dans le milieu indépendamment les unes des autres. En ce qui concerne la propagation dans l'air, on est donc ramené à une théorie scalaire, où l'onde est représentée par une fonction u(P,t). Cette fonction représente en réalité une des six composantes du champ électromagnétique.

En revanche, il n'est pas évident que la théorie scalaire puisse s'appliquer au problème de la diffraction. L'interaction de l'onde avec l'écran, particulièrement au voisinage de l'ouverture, peut en effet introduire un couplage entre les différentes composantes du champ. On se place néanmoins dans le cadre de la théorie scalaire de la diffraction, qui suppose que ces effets de couplage sont négligeables pour le calcul de l'onde diffractée. Cette théorie suppose que l'onde sur la surface de l'ouverture de l'écran (surface plane) est identique à l'onde qui serait présente sur cette surface en l'absence de l'écran et que la fonction d'onde u(P,t) s'annule en dehors de l'ouverture (à l'arrière de l'écran). En d'autres termes, cette théorie ignore l'interaction physique entre l'onde et l'écran lui-même, interaction qui dépend certainement du matériau de l'écran et de la structure microscopique du bord de l'ouverture. Il est généralement admis que cette théorie fournit des résultats conformes aux observations expérimentales lorsque le diamètre de l'ouverture est très grand devant la longueur d'onde. Par exemple, elle peut s'appliquer dans le domaine des ondes visibles pour une ouverture dont la taille est de l'ordre du dixième du millimètre. En revanche, elle ne s'applique pas pour les réseaux de diffraction dont la période est de l'ordre de 10 micromètre ou moins.

Dans le cadre de cette théorie scalaire et puisque l'onde est supposée monochromatique, on utilisera la représentation complexe de u(P,t) :

$$u(P,t)=Re\left(\underline{U}\left(P\right)exp(-j2\pi \nu t)\right)\text{(3)}$$

L'équation des ondes se traduit pour la grandeur complexe par l'équation de Helmholtz :

$$({\nabla }^2+k^2)\underline{U}\left(P\right)=0\text{(4)}$$

où la constante k (le nombre d'onde) est défini par :

$$k=2\pi n\frac{\nu }{c}=\frac{2\pi }{\lambda }\text{(5)}$$

L'intensité de la lumière en un point P est, par définition, la moyenne temporelle de la norme du vecteur de Poynting en ce point. Dans le cadre de la théorie scalaire, l'intensité est proportionnelle à la moyenne temporelle de u²(P,t). Il est d'usage de définir l'intensité par :

$$I\left(P\right)=<u^2(P,t)>\text{(6)}$$

où la notation $<f\left(t\right)>$ désigne la moyenne temporelle sur une période T :

$$<f\left(t\right)>=\frac{1}{T}{\int }_0^{T}f\left(t\right)dt$$

Si l'on adopte la relation $\mathrm{(6)}$, u au carré a les dimensions de l'intensité (énergie par unité de surface). La fonction d'onde u est donc en $W^{\frac{1}{2}}\cdot m^{-1}$. Ce point est cependant sans importance car on s'intéresse aux variations de l'intensité sur un plan et non pas à l'intensité en tant que grandeur énergétique.

Pour une onde monochromatique, on a :

$$I\left(P\right)=\frac{1}{2}\underline{U}(P,t){\underline{U}}^{*}(P,t)=\frac{1}{2}|A\left(P\right){|}^2\text{(7)}$$

2.b. Définition du problème

L'onde monochromatique émise par une source ponctuelle S rencontre un écran plan dans laquelle une ouverture est percée. Le plan Σ situé juste derrière l'écran est muni d'un repère (Oxy), auquel on ajoute un axe (Oz) perpendiculaire au plan. L'objectif est de déterminer $\underline{U}\left(M\right)$ en un point M du plan Π situé à une distance z de l'écran. On note P un point quelconque de l'ouverture, c'est-à-dire un point du plan Σ.

Figure pleine page

Comme mentionné en introduction, on suppose que la largeur de l'ouverture est très grande devant la longueur d'onde (au moins un facteur 100) et que l'écran est assez mince pour que l'interaction de l'onde avec le bord de l'ouverture n'ait qu'un effet négligeable sur l'onde dans le plan Π. En conséquence de ces hypothèses, on admet que la théorie scalaire peut s'appliquer et que l'amplitude complexe $\underline{U}\left(P\right)$ sur le plan Σ au niveau de l'ouverture est exactement celle de l'onde émise par la source qu'il y aurait sur ce plan en l'absence de l'écran. Autrement dit, l'effet physique du bord de l'ouverture est négligé (en réalité seule la comparaison avec l'expérience permet de valider cette hypothèse). Sur la partie du plan Σ située juste derrière l'écran, on suppose que $\underline{U}\left(P\right)=0$. Il s'agit là aussi d'une hypothèse qui peut en réalité être mise en défaut au voisinage du bord de l'ouverture.

2.c. Principe de Huygens-Fresnel

Lorsqu'une onde se propage dans un milieu homogène, l'état de vibration derrière une surface d'onde peut se calculer en considérant des sources fictives secondaires réparties uniformément sur cette surface. Cette propriété, énoncée par Huygens, a été reprise par Fresnel pour traiter le problème de la diffraction. L'hypothèse énoncée par Fresnel est que les points de la surface Σ situés sur l'ouverture sont, du point de vue d'un observateur situé derrière l'écran, équivalents à des sources fictives dont l'état de vibration, c'est-à-dire l'amplitude complexe $\underline{U}\left(P\right)$, est celui de l'onde émise par la source S au point P en l'absence d'écran. Pour un point P de l'ouverture, on a d'après cette hypothèse :

$$\underline{U}\left(P\right)=C\frac{exp\left(jk{r}_s\right)}{r_s}\text{(8)}$$

où r_s désigne la distance de S au point P. La constante C permet de garantir les dimensions correctes pour l'amplitude complexe (C est en $W^{\frac{1}{2}}$).

La démonstration de la validité de l'hypothèse de Fresnel (compte tenu des hypothèses énoncées plus haut) a été faite par Kirchhoff [1]. Une démonstration similaire mais plus rigoureuse due à Rayleigh et Sommerfeld [1] conduit au résultat suivant pour l'amplitude complexe au point M :

$$\underline{U}\left(M\right)=\frac{1}{j\lambda }{\iint }_{\Sigma }\underline{U}\left(P\right)\frac{exp\left(jkr\right)}{r}cos\theta ds\text{(9)}$$

où r désigne la distance entre les points P et M et θ l'angle entre la normale à l'ouverture (le vecteur unitaire ${\overrightarrow{u}}_z$) et le vecteur $\overrightarrow{PM}$.

Ce résultat montre que l'intuition de Fresnel était correcte, mis à part le facteur $\frac{1}{j\lambda }$ et le facteur d'oblicité $cos\theta$. Le résultat obtenu par Kirchoff diffère par le facteur d'oblicité, qui fait intervenir aussi l'angle entre le vecteur $\overrightarrow{SP}$ et la normale.

L'intensité de la lumière au point M est

$$I\left(M\right)=\frac{1}{2}|\underline{U}\left(M\right){|}^2\text{(10)}$$

Il faut remarquer que le facteur 1/j dans l'expression $\mathrm{(9)}$, qui signifie que l'onde émise par P est en quadrature par rapport à l'onde reçue, est sans conséquence sur l'intensité. Fresnel ne pouvait donc pas deviner ce facteur sur la base d'observations. Quant au facteur 1/λ, sa constatation expérimentale est sans doute très délicate et en tout cas impossible à l'époque de Fresnel.

L'intégrale $\mathrm{(9)}$ est généralement considérée comme l'expression mathématique du principe de Huygens-Fresnel, bien que ce dernier n'ait pas prévu le facteur 1/(jλ).

Notons $(\xi ,\eta )$ les coordonnées du point P et (x,y) celles du point M. La distance r s'écrit :

$$r=\sqrt{z^2+{(x-\xi )}^2+{(y-\eta )}^2}\text{(11)}$$

On a par ailleurs :

$$cos\theta =\frac{z}{r}\text{(12)}$$

L'intégrale $\mathrm{(9)}$ s'écrit donc :

$$\underline{U}(x,y)=\frac{z}{j\lambda }{\iint }_{\Sigma }\underline{U}(\xi ,\eta )\frac{exp\left(jkr\right)}{r^2}d\xi d\eta \text{(13)}$$

Rappelons que $\underline{U}(\xi ,\eta )$ est égale à l'amplitude de l'onde émise par la source S et qu'elle est nulle sur la partie de Σ située derrière l'écran.

Introduisons la fonction de transmission de l'écran $\tau (\xi ,\eta )$, vallant 1 sur l'ouverture, 0 sur la partie de Σ située derrière l'écran. On a finalement :

$$\underline{U}(\xi ,\eta )=\tau (\xi ,\eta )C\frac{exp\left(jk{r}_s\right)}{r_s}\text{(14)}$$

où r_s dépend évidemment de $(\xi ,\eta )$ et de la position de S.

2.d. Approximation de Fresnel

La distance r s'écrit :

$$r=z\sqrt{1+(\frac{x-\xi }{z}{)}^2+(\frac{y-\eta }{z}{)}^2}\text{(15)}$$

Les écarts $x-\xi$ et $y-\eta$ sont supposés petits devant z. L'approximation de Fresnel consiste à ne garder que les termes d'ordre 2 :

$$r\approx z(1+\frac{1}{2}(\frac{x-\xi }{z}{)}^2+\frac{1}{2}\left(\frac{y-\eta }{z}{)}^2\right)\text{(16)}$$

L'intégrale $\mathrm{(13)}$ s'écrit :

$$\underline{U}(x,y)=\frac{exp\left(jkz\right)}{j\lambda z}\iint \underline{U}(\xi ,\eta )exp\left[j\frac{k}{2z}\right({(x-\xi )}^2+{(y-\eta )}^2\left)\right]d\xi d\eta \text{(17)}$$

ou encore :

$$\underline{U}(x,y)=\frac{exp\left(jkz\right)}{j\lambda z}exp\left(j\frac{\pi }{\lambda z}(x^2+y^2)\right)\iint \left[\underline{U}(\xi ,\eta )exp\left(j\frac{\pi }{\lambda z}({\xi }^2+{\eta }^2)\right)\right]exp(-j\frac{2\pi }{\lambda z}(x\xi +y\eta ))d\xi d\eta \text{(18)}$$

Afin d'exprimer l'intégrale sous la forme d'une transformée de Fourier, introduisons les fréquences spatiales :

$$\begin{array}{ccc}\hfill & f_x=\frac{x}{\lambda z}\hfill & \hfill \text{(19)}\\ \hfill & f_y=\frac{y}{\lambda z}\hfill & \hfill \text{(20)}\end{array}$$

et la fonction :

$$\underline{V}(\xi ,\eta )=\underline{U}(\xi ,\eta )exp\left(j\frac{\pi }{\lambda z}({\xi }^2+{\eta }^2)\right)\text{(21)}$$

Si l'objectif est de déterminer l'intensité de la lumière sur le plan Π, nous pouvons ignorer le facteur exponentiel en dehors de l'intégrale et écrire :

$$\underline{U}(f_x,f_y)=\frac{1}{\lambda z}\iint \underline{V}(\xi ,\eta )exp(-j2\pi (f_x\xi +f_y\eta ))d\xi d\eta \text{(22)}$$

L'amplitude complexe de l'onde sur le plan Π est donc la transformée de Fourier de la fonction $\underline{V}(\xi ,\eta )$ (multipliée par $\frac{1}{\lambda z}$).

Afin d'expliciter la fonction $\underline{U}(\xi ,\eta )$, notons (x_s,y_s,-z_s) les coordonnées de la source S. On a :

$$r_s=\sqrt{z_s^2+{(x_s-\xi )}^2+{(y_s-\eta )}^2}\text{(23)}$$

et :

$$\underline{U}(\xi ,\eta )=\tau (\xi ,\eta )C\frac{exp\left(jk{r}_s\right)}{r_s}\text{(24)}$$

Un cas particulier intéressant est celui où la source est sur l'axe (Oz), à une distance assez grande pour que $r_s\approx z_s$ (source à l'infini). L'ouverture est alors éclairée par une onde plane en incidence normale. On a dans ce cas :

$$\underline{V}(\xi ,\eta )={\underline{V}}_{\infty }(\xi ,\eta )=C\frac{exp\left(jk{z}_s\right)}{z_s}\tau (\xi ,\eta )exp\left(j\frac{\pi }{\lambda z}({\xi }^2+{\eta }^2)\right)\text{(25)}$$

Dans le cas où la source est à distance finie, on peut utiliser la forme approchée suivante pour la distance r_s qui apparaît dans le facteur $exp\left(jk{r}_s\right)$ :

$$r_s\approx z_s(1+\frac{1}{2}(\frac{x_s-\xi }{z_s}{)}^2+\frac{1}{2}\left(\frac{y_s-\eta }{z_s}{)}^2\right)\text{(26)}$$

et l'approximation $r_s\approx z_s$ pour la distance au dénominateur. On obtient alors :

$$\underline{V}(\xi ,\eta )\approx \frac{exp\left(jk{z}_s(1+\frac{x_s^2}{2}+\frac{ys}{2})\right)}{z_s}\tau (\xi ,\eta )exp\left(j\pi (\frac{1}{z}+\frac{1}{z_s})({\xi }^2+{\eta }^2)\right)exp(-j\frac{2\pi }{\lambda z_s}(x_s\xi +y_s\eta ))\text{(27)}$$

Nous pouvons définir deux fréquences spatiales associées à la position de la source :

$$\begin{array}{ccc}\hfill & s_x=\frac{x_s}{\lambda z_s}\hfill & \hfill \text{(28)}\\ \hfill & s_y=\frac{y_s}{\lambda z_s}\hfill & \hfill \text{(29)}\end{array}$$

L'amplitude de l'onde dans le plan Π s'écrit :

$$\underline{U}(f_x,f_y)=C\frac{exp\left(jk{z}_s(1+\frac{x_s^2}{2}+\frac{ys}{2})\right)}{\lambda z{z}_s}\iint \tau (\xi ,\eta )exp\left(j\frac{\pi }{\lambda }(\frac{1}{z}+\frac{1}{z_s})({\xi }^2+{\eta }^2)\right)exp(-j2\pi ((f_x+s_x)\xi +(f_y+s_y)\eta ))d\xi d\eta \text{(30)}$$

En ignorant le facteur exponentiel en dehors de l'intégrale, qui n'a pas d'effet sur l'intensité, on a finalement :

$$\underline{U}(f_x,f_y)=\frac{C}{\lambda z{z}_s}\iint \tau (\xi ,\eta )exp\left(j\frac{\pi }{\lambda }(\frac{1}{z}+\frac{1}{z_s})({\xi }^2+{\eta }^2)\right)exp(-j2\pi ((f_x+s_x)\xi +(f_y+s_y)\eta ))d\xi d\eta \text{(31)}$$

La fonction dont on calcule la transformée de Fourier est :

$$\underline{V}(\xi ,\tau )=\tau (\xi ,\eta )exp\left(j\frac{\pi }{\lambda }(\frac{1}{z}+\frac{1}{z_s})({\xi }^2+{\eta }^2)\right)\text{(32)}$$

Sa transformée de Fourier est définie par :

$$TF\left(\underline{V}\right)(f_x,f_y)=\iint \underline{V}(\xi ,\tau )exp(-j2\pi (f_x\xi +f_y\eta ))d\xi d\eta \text{(33)}$$

L'amplitude de l'onde sur le plan Π s'en déduit :

$$\underline{U}(f_x,f_y)=\frac{C}{\lambda z{z}_s}TF\left(\underline{V}\right)(f_x+s_x,f_y+s_y)\text{(34)}$$

L'intensité lumineuse est :

$$I(f_x,f_y)=\left(\frac{C}{\lambda z{z}_s}{)}^2\right|TF\left(\underline{V}\right)(f_x+s_x,f_y+s_y){|}^2\text{(35)}$$

Remarquons que $\underline{V}$ est sans dimensions et sa transformée de Fourier a donc les dimensions d'une surface. La constante en facteur permet à l'intensité d'avoir les dimensions d'une puissance par unité de surface. Si on s'intéresse aux variations de l'intensité dans le plan Π, cette constante peut cependant être ignorée.

La figure de diffraction est donnée par la fonction $I(f_x,f_y)$. Voyons l'influence de la position de la source si z_s et fixé. Lorsqu'on change la position de la source, les fréquences s_x et s_y changent mais la fonction $\underline{V}(\xi ,\eta )$ est inchangée. En conséquence, la figure de diffraction pour une position $(s_x,s_y)$ quelconque de la source est identique à la figure de diffraction lorsque la source se trouve sur l'axe mais elle est centrée sur la position définie par f_x=-s_x et f_y=-s_y.

2.e. Approximation de Fraunhofer

Lorsque z et z_s sont très grands devant $k({\xi }^2+{\eta }^2)$, l'approximation de Fraunhofer consiste à négliger la phase dans l'expression $\mathrm{(32)}$ :

$$\underline{V}(\xi ,\eta )=\tau (\xi ,\eta )\text{(36)}$$

On a alors :

$$\underline{U}(f_x,f_y)=\frac{C}{\lambda z{z}_s}\iint \tau (\xi ,\eta )exp(-j2\pi ((f_x+s_x)\xi +(f_y+s_y)\eta ))d\xi d\eta \text{(37)}$$

L'amplitude dans le plan Π est obtenue par transformée de Fourier de la fonction de transmission de l'écran :

$$\underline{U}(f_x,f_y)=\frac{C}{\lambda z{z}_s}TF\left(\tau \right)(f_x+s_x,f_y+s_y)\text{(38)}$$

L'approximation de Fraunhofer donne un résultat exact si S et M sont à l'infini par rapport à l'écran (z_s et z tendent vers l'infini). Expérimentalement, on peut obtenir cette configuration en plaçant S dans le plan focal objet d'une lentille convergente L₁ et le plan Π sur le plan focal image d'une lentille convergente L₂, comme le montre la figure suivante :

Figure pleine page

La projection du vecteur $\overrightarrow{O_{2}M}$ sur le plan (0xy) fait avec l'axe (Oz) un angle α; sa projection sur le plan (Oyz) fait avec l'axe (Oz) un angle β. Ces angles sont reliés aux coordonnées du point M par :

$$\begin{array}{ccc}\hfill & \alpha =\frac{x}{f{\text{'}}_2}\hfill & \hfill \text{(39)}\\ \hfill & \beta =\frac{y}{f{\text{'}}_2}\hfill & \hfill \text{(40)}\end{array}$$

Les fréquences spatiales associées à ce point sont :

$$\begin{array}{ccc}\hfill & f_x=\frac{x}{\lambda f{\text{'}}_2}=\frac{\alpha }{\lambda }\hfill & \hfill \text{(41)}\\ \hfill & f_y=\frac{y}{\lambda f{\text{'}}_2}=\frac{\beta }{\lambda }\hfill & \hfill \text{(42)}\end{array}$$

De manière similaire, la position de la source peut être définie par deux angles :

$$\begin{array}{ccc}\hfill & {\alpha }_s=-\frac{x_s}{f{\text{'}}_1}\hfill & \hfill \text{(43)}\\ \hfill & {\beta }_s=-\frac{y_s}{f{\text{'}}_1}\hfill & \hfill \text{(44)}\end{array}$$

et les fréquences spatiales associées à cette position sont :

$$\begin{array}{ccc}\hfill & s_x=\frac{x_s}{\lambda f{\text{'}}_1}=-\frac{{\alpha }_s}{\lambda }\hfill & \hfill \text{(45)}\\ \hfill & s_y=\frac{y_s}{\lambda f{\text{'}}_1}=-\frac{{\beta }_s}{\lambda }\hfill & \hfill \text{(46)}\end{array}$$

La figure de diffraction pour une position $(s_x,s_y)$ quelconque de la source est identique à la figure de diffraction lorsque la source se trouve sur l'axe mais elle est centrée sur la position définie par f_x=-s_x et f_y=-s_y.

2.f. Transformée de Fourier discrète

Rappelons que l'amplitude diffractée, dans l'approximation de Fresnel ou de Fraunhofer, est proportionnelle à la transformée de Fourier (TF) suivante :

$$\underline{F}(f_x,f_y)=\iint \underline{V}(\xi ,\tau )exp(-j2\pi (f_x\xi +f_y\eta ))d\xi d\eta \text{(47)}$$

La transformée de Fourier discrète (TFD) permet d'obtenir une approximation de la TF. La méthode de calcul est similaire à celle mise en œuvre pour le calcul de la transformée de Fourier d'une image.

Il faut tout d'abord définir sur le plan Σ un domaine rectangulaire de taille L_x,L_y centré en (0,0). Ce domaine doit évidemment contenir l'ouverture toute entière et doit même être beaucoup plus grande que l'ouverture. Plus le domaine est grand, plus la résolution fréquentielle sera grande. On échantillonne le domaine de la manière suivante :

$$\begin{array}{ccc}\hfill & {\xi }_k=-\frac{L_x}{2}+k\frac{L_x}{N_x}0\le k\le N_x-1\hfill & \hfill \text{(48)}\\ \hfill & {\eta }_m=-\frac{L_y}{2}+m\frac{L_y}{N_y}0\le m\le N_y-1\hfill & \hfill \text{(49)}\\ \hfill & {\underline{V}}_{m,k}=\underline{V}({\xi }_m,{\eta }_n)\hfill & \hfill \text{(50)}\end{array}$$

On obtient ainsi une matrice de N_x colonnes et N_y lignes.

Considérons une approximation de la transformée de Fourier obtenue par la méthode des rectangles sur le domaine considéré :

$$\underline{F}(f_x,f_y)\simeq \frac{L_x{L}_y}{N_x{N}_y}exp\left(j\pi (f_x{L}_x+f_y{L}_y)\right)\sum _{k=0}^{N_x-1}\sum _{m=0}^{N_y-1}{\underline{V}}_{m,k}exp(-j2\pi f_{x}k\frac{L_x}{N_x})exp(-j2\pi f_{y}m\frac{L_y}{N_y})\text{(51)}$$

En pratique, on calcule cette approximation pour les fréquences spatiales suivantes :

$$\begin{array}{ccc}\hfill & f_{x,n}=\frac{n}{L_x}0\le n\le N_x-1\hfill & \hfill \text{(52)}\\ \hfill & f_{y,l}=\frac{l}{L_y}0\le l\le N_y-1\hfill & \hfill \text{(53)}\end{array}$$

La transformée de Fourier discrète (TFD) associe à la matrice $\underline{V}$ une matrice $\underline{W}$ de mêmes dimensions, définie par :

$${\underline{W}}_{n,l}=\sum _{k=0}^{N_x-1}\sum _{m=0}^{N_y-1}{\underline{V}}_{m,k}exp(-j2\pi k\frac{n}{N_x})exp(-j2\pi m\frac{l}{N_y})$$

L'approximation de la TF pour les fréquences définies ci-dessus est finalement :

$${\underline{F}}_{n,l}=\frac{L_x{L}_y}{N_x{N}_y}exp\left(j\pi (f_x{L}_x+f_y{L}_y)\right){\underline{W}}_{n,l}\text{(54)}$$

La fréquence d'échantillonnage est constituée par la paire :

$$f_e=(\frac{N_x}{L_x},\frac{N_y}{L_y})$$

On vérifie facilement la propriété suivante de la TFD:

$$\begin{array}{ccc}\hfill & W_{n+N_x,l}=V_{n,l}\hfill & \hfill \text{(55)}\\ \hfill & W_{n,l+N_y}=V_{n,l}\hfill & \hfill \text{(56)}\end{array}$$

Cette propriété correspond à la périodicité du spectre de l'image échantillonnée. La période du spectre est égale à la fréquence d'échantillonnage. En général, on cherche à calculer une approximation de la transformée de Fourier ne faisant pas apparaître cet effet de l'échantillonnage. Pour cela, il faut limiter la plus haute fréquence à la moitié de la fréquence d'échantillonnage. On suppose que N_x et N_y sont des puissances de 2 (une condition nécessaire pour l'utilisation de l'algorithme de transformée de Fourier rapide). La TFD donne une approximation de la transformée de Fourier pour les fréquences suivantes :

$$\begin{array}{ccc}\hfill & f_{x,n}=\frac{n}{L_x}0\le n\le \frac{N_x}{2}\hfill & \hfill \text{(57)}\\ \hfill & f_{y,l}=\frac{l}{L_y}0\le l\le \frac{N_y}{2}\hfill & \hfill \text{(58)}\end{array}$$

De plus, la périodicité de la TFD permet d'accéder aux opposées de ces fréquences. En effet :

$$\begin{array}{ccc}\hfill & V_{-n,l}=V_{-n+N_x,l}\hfill & \hfill \text{(59)}\\ \hfill & V_{n,-l}=V_{n,-l+N_y}\hfill & \hfill \text{(60)}\end{array}$$

Les N_x indices n correspondent aux fréquences suivantes :

$$0,\frac{1}{L_x},\frac{2}{L_x},\dots ,\frac{\frac{N_x}{2}}{L_x},-\frac{\frac{N_x}{2}-1}{L_x},\dots ,-\frac{2}{L_x},-\frac{1}{L_x}$$

La même relation est valable pour les fréquences de l'axe y.

La fréquence nulle correspond à un point M sur l'axe optique. On doit placer ce point au milieu de la matrice, de manière à obtenir des fréquences croissantes en fonction des indices. Dans cette matrice, le point de fréquence nulle doit être placé aux indices suivants :

$$(\frac{N_y}{2}-1,\frac{N_x}{2}-1)$$

Pour l'axe x, il y a donc $\frac{N_x}{2}+1$ indices correspondant à des fréquences positives (ou nulle) et $\frac{N_x}{2}-1$ indices correspondant à des fréquences strictement négatives.

Pour construire la matrice où la fréquence nulle est centrée, il faut utiliser la transformation d'indice suivante :

$$\begin{array}{cc}\hfill & n_c=\frac{N_x}{2}-1+n\text{si}0\le n\le \frac{N_x}{2}\hfill \\ \hfill & n_c=-\frac{N_x}{2}-1+n\text{si}\frac{N_x}{2}+1\le n\le N_x-1\hfill \end{array}$$

Considérons par exemple le calcul de la figure de diffraction d'une ouverture circulaire dans l'approximation de Fraunhofer.

import numpy as np
from matplotlib.pyplot import *
R = 0.1 # rayon de l'ouverture en mm
def tau(xi,eta):
    if xi**2+eta**2 < R**2: return 1.0
    else: return 0.0
tau_ufunc = np.frompyfunc(tau,2,1)
Lx = 100*R
Ly = 100*R 
Nx = 2**11
Ny = 2**11
xi,eta  = np.meshgrid(-Lx/2+np.arange(Nx)*Lx/Nx,-Ly/2+np.arange(Ny)*Ly/Ny)
V = tau_ufunc(xi,eta)
V = V.astype(float)
figure()
imshow(V,vmin=0,vmax=1.0,cmap='gray',extent=[-Lx,Lx,-Ly,Ly])
xlabel(r"$\xi\ (\rm mm)$",fontsize=16)
ylabel(r"$\eta\ (\rm mm)$",fontsize=16)
            

fig1.pdf

Le calcul de la TFD se fait avec la fonction fft2. La fonction shiftfft permet de placer la fréquence nulle au centre.

from numpy.fft import fft2, fftshift
F = fftshift(fft2(V))*Lx*Lx/(Nx*Ny)
I = np.power(np.absolute(F),2)
figure()
fx_max = Nx/(2*Lx)
fy_max = Ny/(2*Ly)
imshow(I,cmap='gray',extent=[-fx_max,fx_max,-fy_max,fy_max])
xlabel(r"$f_x\ (\rm mm^{-1})$",fontsize=16)
ylabel(r"$f_y\ (\rm mm^{-1})$",fontsize=16)
            

fig2.pdf

Avec l'échelle complète de fréquences, la tâche de diffraction est très petite. Voici une vue rapprochée :

xlim(-20,20)
ylim(-20,20)
            

fig3.pdf

Voici le profil d'intensité sur une ligne passant par le centre :

Ix = I[Ny//2,:]
fx = -fx_max + np.arange(Nx)*2*fx_max/Nx
figure()
plot(fx,Ix)
grid()
xlim(-20,20)
xlabel(r"$f_x\ (\rm mm^{-1})$",fontsize=16)
ylabel(r"$I$",fontsize=16)
            

fig4.pdf

Pour obtenir un angle en radians en échelle horizontale, il suffit de multiplier la fréquence par la longueur d'onde :

lamb = 500e-6 # longueur d'onde en mm
alpha = fx*lamb
figure()
plot(alpha*1e3,Ix)
grid()
xlim(-20*lamb*1e3,20*lamb*1e3)
xlabel(r"$\alpha\ (\rm mrad)$",fontsize=16)
ylabel(r"$I$",fontsize=16)
            

fig5.pdf

Voyons un deuxième exemple : la diffraction par deux trous circulaires (expérience de Young) :

import numpy as np
from matplotlib.pyplot import *
R = 0.1 # rayon d'un trou en mm
a = 0.5 # distance entre les trous
def tau(xi,eta):
    if (xi-a/2)**2+eta**2 < R**2 or (xi+a/2)**2+eta**2 < R**2 : return 1.0
    else: return 0.0
tau_ufunc = np.frompyfunc(tau,2,1)
Lx = 100*R
Ly = 100*R 
Nx = 2**11
Ny = 2**11
xi,eta  = np.meshgrid(-Lx/2+np.arange(Nx)*Lx/Nx,-Ly/2+np.arange(Ny)*Ly/Ny)
V = tau_ufunc(xi,eta)
V = V.astype(float)
figure()
imshow(V,vmin=0,vmax=1.0,cmap='gray',extent=[-Lx,Lx,-Ly,Ly])
xlabel(r"$\xi\ (\rm mm)$",fontsize=16)
ylabel(r"$\eta\ (\rm mm)$",fontsize=16)
            

fig6.pdf

F = fftshift(fft2(V))*Lx*Lx/(Nx*Ny)
I = np.power(np.absolute(F),2)
figure()
fx_max = Nx/(2*Lx)
fy_max = Ny/(2*Ly)
imshow(I,cmap='gray',extent=[-fx_max,fx_max,-fy_max,fy_max])
xlabel(r"$f_x\ (\rm mm^{-1})$",fontsize=16)
ylabel(r"$f_y\ (\rm mm^{-1})$",fontsize=16)
xlim(-20,20)
ylim(-20,20)
            

fig7.pdf

Ix = I[Ny//2,:]
fx = -fx_max + np.arange(Nx)*2*fx_max/Nx
figure()
plot(fx,Ix)
grid()
xlim(-20,20)
xlabel(r"$f_x\ (\rm mm^{-1})$",fontsize=16)
ylabel(r"$I$",fontsize=16)
            

fig8.pdf

Voyons à présent comment calculer la figure de diffraction pour une position quelconque de la source, définie par les fréquences spatiales s_x,s_y. La nouvelle figure de diffraction est centrée en f_x=-s_x,f_y=-s_y. Il faut donc appliquer un décalage sur les éléments de la matrice car, pour une raison qui apparaîtra plus loin, la taille de la matrice ne doit pas changer. Les nombres de pixels de décalage à appliquer horizontalement et verticalement sont :

$$\begin{array}{ccc}\hfill & \Delta i=-E\left[s_x{L}_x\right]\hfill & \hfill \text{(61)}\\ \hfill & \Delta j=-E\left[s_y{L}_y\right]\hfill & \hfill \text{(62)}\end{array}$$

Bien évidemment, le décalage ne doit pas être trop grand afin qu'aucun pixel dont l'intensité est significative ne soit perdu dans l'opération. Il faut faire ce décalage le plus efficacement possible (en terme de vitesse) et donc éviter d'écrire une boucle pour décaler pixels par pixels. La fonction np.roll permet d'effectuer une rotation des éléments. Dans le cas d'une rotation d'une unité vers la gauche, la première rangée verticale de pixels se retrouve tout à droite. Cela ne devrait pas poser problème dans la mesure où les pixels proches des bords de l'image ont une intensité pratiquement nulle (la fréquence d'échantillonnage doit être assez grande pour cela). La fonction suivante effectue le décalage :

def decalage(I,sx,sy,Lx,Ly):
    (Ny,Nx) = I.shape
    Di = -int(sx*Lx)
    Dj = -int(sy*Ly)
    I1 = np.roll(I,Dj,axis=0)
    return np.roll(I1,Di,axis=1)
                

Voici un décalage appliqué à la figure de diffraction précédente :

I1 = decalage(I,5,0,Lx,Ly)
figure()
fx_max = Nx/(2*Lx)
fy_max = Ny/(2*Ly)
imshow(I1,cmap='gray',extent=[-fx_max,fx_max,-fy_max,fy_max])
xlabel(r"$f_x\ (\rm mm^{-1})$",fontsize=16)
ylabel(r"$f_y\ (\rm mm^{-1})$",fontsize=16)
xlim(-20,20)
ylim(-20,20)
                

fig9.pdf

3.a. Lumière quasi monochromatique

Lorsque la source de lumière est étendue spatialement, le caractère non parfaitement monochromatique de la lumière ne peut être ignoré. En effet, deux points distincts de la source étendue émettent des vibrations incohérentes entre elles.

Lorsque la lumière est polychromatique, une composante du champ électromagnétique en un point P s'écrit [2] :

$$u(P,t)={\int }_{-\infty }^{\infty }\tilde{\underline{U}}(P,\nu )exp(-j2\pi \nu t)d\nu \text{(63)}$$

La fonction $\tilde{\underline{U}}(P,\nu )$ définit le spectre en fréquence de la lumière. Puisque $u(P,t)$ est réel, on a :

$$\underline{\tilde{U}}(P,-\nu )={\underline{\tilde{U}}}^{*}(P,\nu )\text{(64)}$$

$u(P,t)$ est la partie réelle d'une fonction à valeurs complexes (appelée signal analytique associé à u(P,t)) :

$$\underline{U}(P,t)=2{\int }_0^{\infty }\tilde{\underline{U}}(P,\nu )exp(-j2\pi \nu t)d\nu \text{(65)}$$

L'intensité lumineuse est définie par :

$$I\left(P\right)=\underset{T\to \infty }{lim}\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}{u}^2(P,t)dt\text{(66)}$$

L'intensité s'exprime en fonction de l'amplitude complexe par la relation [2] :

$$I\left(P\right)=\frac{1}{2}<|\underline{U}(P,t){|}^2>=\frac{1}{2}<\underline{U}(P,t){\underline{u}}^{*}(P,t)>\text{(67)}$$

où la notation $<f\left(t\right)>$ désigne la moyenne temporelle de la fonction f(t), c'est-à-dire :

$$<f\left(t\right)>=\underset{T\to \infty }{lim}\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f\left(t\right)dt\text{(68)}$$

Remarquons que, contrairement au cas d'une onde monochromatique, la fonction $|\underline{U}(P,t){|}^2$ est en général non périodique et dépend du temps.

Par définition, le spectre d'une lumière quasi monochromatique, c'est-à-dire $\tilde{\underline{U}}(P,\nu )$, ne prend des valeurs non nulles que pour une petite bande de largeur $\Delta \nu$ centrée sur $\overline{\nu }$ (on a $\Delta \nu \ll \overline{\nu }$).

La fonction d'auto-corrélation du champ en un point P est définie par :

$$\Gamma (P,P,\tau )=<\underline{U}(P,t+\tau ){\underline{U}}^{*}(P,t)>\text{(69)}$$

La fonction d'auto-corrélation contient toute l'information utile sur la cohérence temporelle. Elle s'annule pratiquement lorsque $\tau$ est supérieur au temps de cohérence t_c. Le temps de cohérence est environ égal à l'inverse de la largeur de bande de fréquence :

$$t_c\approx \frac{1}{\Delta \nu }\text{(70)}$$

Pour deux points P₁ et P₂ distincts, par exemple deux points distincts de l'ouverture, on définit la fonction de corrélation suivante :

$$\Gamma (P_1,P_2,\tau )=<\underline{U}(P_1,t+\tau ){\underline{U}}^{*}(P_2,t)>\text{(71)}$$

Cette fonction de corrélation est aussi nommée fonction de cohérence mutuelle. Elle contient toute l'information utile sur la cohérence temporelle mais aussi sur la cohérence spatiale (cohérence entre deux points distincts pour $\tau =0$).

En ce qui concerne la source étendue, on considère que deux points distincts S₁ et S₂ de la source émettent de manière parfaitement incohérente, c'est-à-dire :

$$\Gamma (S_1,S_2,\tau )=0\forall \tau \text{(72)}$$

Il est intéressant de normaliser la fonction de corrélation entre deux points :

$$\gamma (P_1,P_2,\tau )=\frac{\Gamma (P_1,P_2,\tau )}{\sqrt{I_1{I}_2}}\text{(73)}$$

avec :

$$\begin{array}{ccc}\hfill & I_1=\frac{1}{2}<\underline{U}(P_1,t){\underline{U}}^{*}(P_1,t)>\hfill & \hfill \text{(74)}\\ \hfill & I_2=\frac{1}{2}<\underline{U}(P_2,t){\underline{U}}^{*}(P_2,t)>\hfill & \hfill \text{(75)}\end{array}$$

Cette fonction prend des valeurs dont le module est dans l'intervalle [0,1]. Si $\left|\gamma (P_1,P_2,0)\right|=1$, les ondes aux points P₁ et P₂ sont dites cohérentes. Si $\left|\gamma (P_1,P_2,0)\right|=0$, les ondes aux points P₁ et P₂ sont incohérentes. Pour une valeur intermédiaire entre 0 et 1, on parle de cohérence partielle. Qu'il y a ait cohérence parfaite ou cohérence partielle, le module de la fonction de corrélation décroît avec $\tau$ et s'annule au delà du temps de cohérence.

3.b. Théorème de Van Cittert-Zernike

La théorie de la diffraction exposée précédemment s'applique si l'ouverture est éclairée de manière cohérente ou quasi cohérente, c'est-à-dire si pour deux points quelconques P₁ et P₂ sur l'ouverture on a $\left|\gamma (P_1,P_2,0)\right|\approx 1$.

Le théorème de Van Cittert-Zernike [2] permet de déterminer la fonction de corrélation $\gamma (P_1,P_2,0)$ lorsque la source est étendue. Ce théorème repose sur l'hypothèse (déjà donnée plus haut) que deux points distincts de la source émettent de manière incohérente :

$$\Gamma (S_1,S_2,0)=0\text{(76)}$$

Notons S un point quelconque de la source étendue et σ la surface de cette source (supposée parallèle à l'ouverture). Définissons les distances R₁=SP₁ et R₂=SP₂. Soit I(S) l'intensité par unité d'aire au points S.

Figure pleine page

Le théorème de Van Cittert-Zernike établit que :

$$\gamma (P_1,P_2,0)=\frac{1}{\sqrt{I_1{I}_2}}{\iint }_{\sigma }I\left(S\right)\frac{exp\left(j\overline{k}(R_1-R_2)\right)}{R_1{R}_2}dS\text{(77)}$$

où $\overline{k}$ désigne le nombre d'onde moyen :

$$\overline{k}=2\pi n\frac{\overline{\nu }}{c}=\frac{2\pi }{\overline{\lambda }}\text{(78)}$$

et :

$$\begin{array}{ccc}\hfill & I_1={\int }_{\sigma }\frac{I\left(S\right)}{R_1^2}dS\hfill & \hfill \text{(79)}\\ \hfill & I_2={\int }_{\sigma }\frac{I\left(S\right)}{R_2^2}dS\hfill & \hfill \text{(80)}\end{array}$$

Si on considère le cas d'une source d'émission uniforme, on a :

$$\gamma (P_1,P_2,0)=\frac{I\left(S\right)}{\sqrt{I_1{I}_2}}{\iint }_{\sigma }\frac{exp\left(j\overline{k}(R_1-R_2)\right)}{R_1{R}_2}dS\text{(81)}$$

Le calcul de cette intégrale pour une source de surface σ donnée est similaire au calcul de l'amplitude de l'onde diffractée par une ouverture (exposé précédemment). La surface σ est similaire à la surface de l'ouverture et les points P₁ et P₂ sont similaires aux points S et M considérés dans le calcul de la diffraction avec une source ponctuelle. La seule différence est qu'il n'y a pas de facteur d'oblicité.

Explicitons les distances en fonction des coordonnées (x_s,y_s,-z_s) du points de la source étendu et des coordonnées des points P₁ et P₂ situés sur l'ouverture :

$$\begin{array}{ccc}\hfill & R_1=\sqrt{z_s^2+{(x_s-{\xi }_1)}^2+{(y_s-{\eta }_1)}^2}\hfill & \hfill \text{(82)}\\ \hfill & R_2=\sqrt{z_s^2+{(x_s-{\xi }_2)}^2+{(y_s-{\eta }_2)}^2}\hfill & \hfill \text{(83)}\end{array}$$

Utilisons pour ces deux distances la forme approchée suivante, valable lorsque la source est à une distance grande devant sa propre taille et devant celle de l'ouverture :

$$\begin{array}{ccc}\hfill & R_1\approx z_s+\frac{{(x_s-{\xi }_1)}^2+{(y_s-{\eta }_1)}^2}{2z_s}\hfill & \hfill \text{(84)}\\ \hfill & R_2\approx z_s+\frac{{(x_s-{\xi }_2)}^2+{(y_s-{\eta }_2)}^2}{2z_s}\hfill & \hfill \text{(85)}\end{array}$$

On en déduit :

$$R_1-R_2=\frac{{\xi }_1^2+{\eta }_1^2-{\xi }_2^2-{\eta }_2^2}{2z_s}+\frac{x_s({\xi }_2-{\xi }_1)+y_s({\eta }_2-{\eta }_1)}{z_s}\text{(86)}$$

Posons :

$$\begin{array}{ccc}\hfill & p=\frac{{\xi }_1-{\xi }_2}{z_s}\hfill & \hfill \text{(87)}\\ \hfill & q=\frac{{\eta }_1-{\eta }_2}{z_s}\hfill & \hfill \text{(88)}\\ \hfill & \psi =\frac{\overline{k}({\xi }_1^2+{\eta }_1^2-{\xi }_2^2-{\eta }_2^2)}{2z_s}\hfill & \hfill \text{(89)}\end{array}$$

Dans l'intégrale $\mathrm{(81)}$, on applique l'approximation $R_1{R}_2\approx z_s^2$. Si σ désigne aussi l'aire de la surface de la source, on obtient :

$$\gamma (P_1,P_2,0)=\frac{exp\left(j\psi \right)}{\sigma }{\iint }_{\sigma }exp(-j\overline{k}(p{x}_s+q{y}_s))d{x}_{s}d{y}_s\text{(90)}$$

Introduisons les fréquence spatiales suivantes :

$$\begin{array}{ccc}\hfill & f_p=\frac{p}{\lambda }=\frac{{\xi }_1-{\xi }_2}{\lambda z_s}\hfill & \hfill \text{(91)}\\ \hfill & f_q=\frac{p}{\lambda }=\frac{{\eta }_1-{\eta }_2}{\lambda z_s}\hfill & \hfill \text{(92)}\end{array}$$

Il vient :

$$\gamma (P_1,P_2,0)=\frac{exp\left(j\psi \right)}{\sigma }{\iint }_{\sigma }exp(-j2\pi (f_p{x}_s+f_q{y}_s))d{x}_{s}d{y}_s\text{(93)}\text{(94)}$$

La fonction de corrélation entre les ondes aux points P₁ et P₂ est finalement la transformée de Fourier de la fonction de transmission définissant la source. Il est possible de calculer précisément cette fonction de corrélation pour des formes de source simples (disque ou rectangle). Nous pouvons aussi utiliser des propriétés générales de la transformée de Fourier. Si L désigne l'ordre de grandeur de la largeur de la source, $\left|\gamma (P_1,P_2,0)\right|$ est proche de 1 si les fréquences spatiales f_p et f_q sont très petites devant 1/L. On en conclut que les deux points P₁ et P₁ sont cohérents si :

$$\begin{array}{ccc}\hfill & |{\xi }_1-{\xi }_2|\ll \frac{\overline{\lambda }{z}_s}{L}\hfill & \hfill \text{(95)}\\ \hfill & |{\eta }_1-{\eta }_2|\ll \frac{\overline{\lambda }{z}_s}{L}\hfill & \hfill \text{(96)}\end{array}$$

Si ces conditions sont vérifiées pour toutes paires de points de l'ouverture, les différents points de l'ouverture peuvent être considérés comme cohérents (éclairage de l'ouverture cohérent) et la théorie développée précédemment dans l'hypothèse d'une source ponctuelle monochromatique donne un résultat correct en ce qui concerne la figure de diffraction. Si au contraire une de ces deux conditions n'est pas vérifiée, différents points de l'ouverture peuvent présenter une cohérence partielle ou même pas du tout de cohérence s'ils sont très éloignés. On peut dire dans ce cas que l'éclairage de l'ouverture n'est pas cohérent.

On peut aussi voir les deux inégalités précédentes comme des conditions pour que la source puisse être considérée comme ponctuelle. En effet, une source peut être considérée comme ponctuelle si elle produit un éclairage cohérent sur l'ouverture.

Donnons l'expression exacte de $\gamma (P_1,P_2,0)$ dans le cas d'une source de forme rectangulaire. Soient l_x,l_y les dimensions de la source. Le calcul de la transformée de Fourier d'une fonction de transmission rectangulaire présenté en annexe (diffraction par une ouverture rectangulaire) permet d'établir que :

$$\gamma (P_1,P_2,0)=exp\left(j\psi \right)sinc\left(\pi f_p{l}_x\right)sinc\left(\pi f_q{l}_y\right)\text{(97)}$$

Explicitement en fonction des coordonnées des points P₁ et P₂ :

$$\gamma (P_1,P_2,0)=exp\left(j\frac{\pi }{\overline{\lambda }}\frac{({\xi }_1^2+{\eta }_1^2-{\xi }_2^2-{\eta }_2^2)}{z_s}\right)sinc\left(\pi \frac{({\xi }_1-{\xi }_2)l_x}{\overline{\lambda }{z}_s}\right)sinc\left(\pi \frac{({\eta }_1-{\eta }_2)l_y}{\overline{\lambda }{z}_s}\right)\text{(98)}$$

Le module de cette fonction d'autocorrélation est proche de 1 pour toutes les paires de points de l'ouverture si les deux conditions suivantes sont vérifiées :

$$\begin{array}{ccc}\hfill & |{\xi }_1-{\xi }_2{|}_{max}\ll \frac{\overline{\lambda }{z}_s}{l_x}\hfill & \hfill \text{(99)}\\ \hfill & |{\eta }_1-{\eta }_2{|}_{max}\ll \frac{\overline{\lambda }{z}_s}{l_y}\hfill & \hfill \text{(100)}\end{array}$$

ce qui nous donne la condition de cohérence spatiale pour l'éclairage d'une ouverture par une source de forme rectangulaire.

3.c. Diffraction en éclairage incohérent

Lorsque la source est très étendue, ou lorsque l'ouverture est très grande, les points de l'ouverture ne sont pas tous cohérents entre eux.

D'après l'expression intégrale $\mathrm{(9)}$ de l'amplitude, l'intensité lumineuse au point M s'écrit :

$$I\left(M\right)=\frac{1}{2{\lambda }^2}<\left({\iint }_{\Sigma }\underline{U}(P,t)\frac{exp\left(j\overline{k}r\right)}{r}cos\theta ds\right)\left({\iint }_{\Sigma }{\underline{U}}^{*}(P\text{'},t)\frac{exp(-j\overline{k}r\text{'})}{r\text{'}}cos\theta \text{'}ds\text{'}\right)>\text{(101)}$$

Ce produit de deux intégrales double peut aussi s'écrire comme une intégrale quadruple :

$$I\left(M\right)=\frac{1}{2{\lambda }^2}{\iint }_{\Sigma }ds{\iint }_{\Sigma }ds\text{'}<\underline{U}(P,t){\underline{U}}^{*}(P\text{'},t)>\frac{exp\left(j\overline{k}(r-r\text{'})\right)}{rr\text{'}}cos\theta cos\theta \text{'}\text{(102)}$$

Cette expression fait intervenir la fonction de corrélation définie précédemment :

$$I\left(M\right)=\frac{1}{2{\lambda }^2}{\iint }_{\Sigma }ds{\iint }_{\Sigma }ds\text{'}\Gamma (P,P\text{'},0)\frac{exp\left(j\overline{k}(r-r\text{'})\right)}{rr\text{'}}cos\theta cos\theta \text{'}\text{(103)}$$

Comme expliqué plus haut, la fonction de corrélation $\Gamma (P,P\text{'},0)$ peut être calculée par transformée de Fourier de la fonction de transmission de la source. Pour une source de forme simple (rectangulaire au circulaire), une expression analytique peut être obtenue. Par exemple, $\mathrm{(98)}$ est l'expression de la fonction de corrélation normalisée pour une source de forme rectangulaire.

Dans le cadre de l'approximation de Fraunhofer et dans le cas où s_x=s_y=0, l'intégrale s'écrit :

$$I\left(M\right)=\frac{C{I}_0}{\lambda z{z}_s}\int d{\xi }_1\int d{\xi }_2\int d{\eta }_1\int d{\eta }_2\tau ({\xi }_1,{\eta }_1)\tau ({\xi }_2,{\eta }_2)\gamma (P_1,P_2,0)exp(-j2\pi (f_x({\xi }_1-{\xi }_2)+f_y({\eta }_1-{\eta }_2)))\text{(104)}$$

où I₀ est une constante homogène à une intensité. Dans le cas d'une ouverture de forme rectangulaire et d'une source de forme rectangulaire, les variables $\xi$ et $\eta$ sont séparables. On obtient le produit de deux intégrales similaires. Par exemple, l'intégrale sur les variables ${\xi }_1,{\xi }_2$ s'écrit :

$${\int }_{-\frac{a}{2}}^{\frac{a}{2}}d{\xi }_1{\int }_{-\frac{a}{2}}^{\frac{a}{2}}d{\xi }_{2}exp\left(j\frac{\pi }{\overline{\lambda }}\frac{{\xi }_1^2-{\xi }_2^2}{z_s}\right)sinc\left(\pi \frac{({\xi }_1-{\xi }_2)L_x}{\overline{\lambda }{z}_s}\right)exp(-j2\pi f_x({\xi }_1-{\xi }_2))\text{(105)}$$

Même dans ce cas simple, les variables ${\xi }_1$ et ${\xi }_2$ ne sont pas séparables (elles ne le sont pas dans le sinus cardinal) et nous ne connaisons pas d'expression analytique de cette intégrale.

L'expression $\mathrm{(103)}$ montre la relation entre la fonction de corrélation et l'intensité mais elle est de peu d'intérêt pour le calcul numérique de l'intensité dans le cas général (et même dans le cas simple évoqué ci-dessus).

Une méthode plus simple et plus générale consiste à utiliser l'hypothèse fondamentale énoncée pour établir le théorème de Van Cittert-Zernike : deux points de la source étendue sont incohérents. Notons $I(s_x,s_y,f_x,f_y)$ l'intensité obtenue avec une source ponctuelle, pour une position quelconque de celle-ci. On calcule tout d'abord l'intensité lorsque s_x=s_y=0.

$$I(0,0,f_x,f_y)=\left(\frac{C}{\lambda z{z}_s}{)}^2\right|TF\left(\underline{V}\right)(f_x,f_y){|}^2\text{(106)}$$

L'intensité pour une position quelconque de la source s'en déduit (à condition de ne pas changer la distance du plan de la source) :

$$I(s_x,s_y,f_x,f_y)=I{(0,0,f_x+s_x,f_y+s_y)}_{}\text{(107)}$$

Rappelons que dans le cadre de l'approximation de Fraunhofer, $\underline{V}$ est simplement la fonction de transmission de l'écran, égale à 1 sur l'ouverture, à 0 sur la partie opaque.

Puisque les différents points de la source sont incohérents, l'intensité avec la source étendue est l'intégrale de cette intensité sur la surface de la source :

$$I(f_x,f_y)=\frac{1}{\sigma }{\iint }_{\sigma }I(s_x,s_y,f_x,f_y)d{s}_{x}d{s}_y\text{(108)}$$

Même dans le cas le plus simple (source et ouverture rectangulaire), il est très difficile (cela nous semble impossible) de calculer cette intégrale analytiquement. Nous pouvons en revanche la calculer numériquement, ce qui permet évidemment de traiter le cas général. Il faut remarquer que le calcul numérique de l'intensité se fait très efficacement par l'agorithme de transformée de Fourier rapide. Il suffit de faire le calcul de $I(0,0,f_x,f_y)$ puisque l'intensité pour une position quelconque de la source s'en déduit par une simple translation de $(s_x,s_y)$. Le calcul numérique de l'intégrale sur la surface de la source $\mathrm{(108)}$ se fera, dans le cas général, par une méthode de Monte-Carlo. Dans le cas d'une intégrale double, une méthode de Monte-Carlo n'apporte pas de gain d'efficacité par rapport à une méthode de quadrature (par exemple la méthode des rectangles). Cependant, cette dernière nécessite de réaliser un quadrillage de la surface de la source, ce qui peut être difficile pour des formes complexes (même pour un disque, la tâche n'est pas simple). Par ailleurs, la méthode de Monte-Carlo a l'avantage de fournir rapidement une approximation grossière de l'intégrale, qui s'affine progressivement au cours du temps de calcul (lorsque le nombre de tirages augmente). La complexité globale du calcul est celle d'une FFT à deux dimensions suivie d'une intégrale double. Pour comparaison, l'application de la relation $\mathrm{(103)}$ nécessiterait tout d'abord le calcul de la fonction de corrélation par une FFT à deux dimensions, suivie d'une intégrale quadruple.

3.d. Calcul numérique

Reprenons l'exemple de la diffraction par deux trous circulaires (expérience de Young). Nous avons une matrice contenant la figure de diffraction pour une source centrée et une fonction de décalage, permettant d'obtenir la figure de diffraction pour une position quelconque de la source. La source est un disque de rayon R_s situé à une distance z_s de l'écran contenant les deux trous. Pour un point de coordonnées x_s,y_s sur le disque de la source, les fréquences spatiales sont :

$$\begin{array}{ccc}\hfill & s_x=\frac{x_s}{\lambda z_s}\hfill & \hfill \text{(109)}\\ \hfill & s_y=\frac{y_s}{\lambda z_s}\hfill & \hfill \text{(110)}\end{array}$$

Le tirage d'un point aléatoire sur le disque de la source, avec une densité de probabilité uniforme, se fait sans difficulté en coordonnées polaires :

$$\begin{array}{cc}\hfill & \theta =2\pi u_1\hfill \\ \hfill & r=R_s\sqrt{u_2}\hfill \end{array}$$

où u₁ et u₁ sont deux nombres aléatoires de densité de probabilité uniforme sur l'intervalle [0,1].

La fonction suivante effectue un tirage aléatoire d'un point sur le disque et renvoie les deux fréquences spatiales :

import numpy.random as random
def tirageSource(Rs,zs,lamb):
    u1 = random.uniform()
    u2 = random.uniform()
    theta = 2*np.pi*u1
    r = Rs*np.sqrt(u2)
    xs = r*np.cos(theta)
    ys = r*np.sin(theta)
    return (xs/(zs*lamb),ys/(zs*lamb))
                   

Pour chaque tirage, on calcule la figure de diffraction correspondant à cette position de la source ponctuelle puis on ajoute l'intensité à l'intensité déjà calculée :

Nt = 100 # nombre de tirage
Rs = 0.02 # en mm
zs = 100
lamb = 500e-6
I2 = I.copy()
for t in range(Nt):
    (sx,sy) = tirageSource(Rs,zs,lamb)
    I2 += decalage(I,sx,sy,Lx,Ly)
figure()
imshow(I2,cmap='gray',extent=[-fx_max,fx_max,-fy_max,fy_max])
xlabel(r"$f_x\ (\rm mm^{-1})$",fontsize=16)
ylabel(r"$f_y\ (\rm mm^{-1})$",fontsize=16)
xlim(-20,20)
ylim(-20,20)
                

fig10.pdf

Ix = I2[Ny//2,:]
fx = -fx_max + np.arange(Nx)*2*fx_max/Nx
figure()
plot(fx,Ix)
grid()
xlim(-20,20)
xlabel(r"$f_x\ (\rm mm^{-1})$",fontsize=16)
ylabel(r"$I$",fontsize=16)
            

fig11.pdf

Pour cette taille de la source, la diminution du contraste des franges est bien visible sur le profil d'intensité. Voici le résultat pour une source plus grande :

Nt = 200 # nombre de tirage
Rs = 0.04 # en mm
zs = 100
lamb = 500e-6
I2 = I.copy()
for t in range(Nt):
    (sx,sy) = tirageSource(Rs,zs,lamb)
    I2 += decalage(I,sx,sy,Lx,Ly)
figure()
imshow(I2,cmap='gray',extent=[-fx_max,fx_max,-fy_max,fy_max])
xlabel(r"$f_x\ (\rm mm^{-1})$",fontsize=16)
ylabel(r"$f_y\ (\rm mm^{-1})$",fontsize=16)
xlim(-20,20)
ylim(-20,20)
                

fig12.pdf

Ix = I2[Ny//2,:]
fx = -fx_max + np.arange(Nx)*2*fx_max/Nx
figure()
plot(fx,Ix)
grid()
xlim(-20,20)
xlabel(r"$f_x\ (\rm mm^{-1})$",fontsize=16)
ylabel(r"$I$",fontsize=16)
            

fig13.pdf

puis pour une source beaucoup plus grande, pour laquelle les deux trous sont éclairés de manière incohérentes :

Nt = 500 # nombre de tirage
Rs = 0.1 # en mm
zs = 100
lamb = 500e-6
I2 = I.copy()
for t in range(Nt):
    (sx,sy) = tirageSource(Rs,zs,lamb)
    I2 += decalage(I,sx,sy,Lx,Ly)
figure()
imshow(I2,cmap='gray',extent=[-fx_max,fx_max,-fy_max,fy_max])
xlabel(r"$f_x\ (\rm mm^{-1})$",fontsize=16)
ylabel(r"$f_y\ (\rm mm^{-1})$",fontsize=16)
xlim(-20,20)
ylim(-20,20)
                

fig14.pdf

Ix = I2[Ny//2,:]
fx = -fx_max + np.arange(Nx)*2*fx_max/Nx
figure()
plot(fx,Ix)
grid()
xlim(-20,20)
xlabel(r"$f_x\ (\rm mm^{-1})$",fontsize=16)
ylabel(r"$I$",fontsize=16)
            

fig15.pdf

Lorsque la source est encore plus grande, le nombre de tirage doit être augmenté.

4.a. Ouverture rectangulaire

Cette partie montre le calcul de l'amplitude diffractée dans le cas d'une ouverture rectangulaire et dans le cadre de l'approximation de Fraunhofer. L'ouverture est un rectangle de largeur a et de longueur b. La fonction de transmission est définie par :

$$\begin{array}{cc}\hfill & \tau (\xi ,\eta )=1si-\frac{a}{2}<\xi <\frac{a}{2}et-\frac{b}{2}<\eta <\frac{b}{2}\hfill \\ \hfill & \tau (\xi ,\eta )=0sinon\hfill \end{array}$$

L'intégrale $\mathrm{(37)}$ s'écrit, dans le cas s_x=s_y=0 :

$$\underline{U}(f_x,f_y)=\frac{C}{\lambda z{z}_s}{\int }_{-\frac{a}{2}}^{\frac{a}{2}}{\int }_{-\frac{b}{2}}^{\frac{b}{2}}exp(-j2\pi (f_x\xi +f_y\eta ))d\xi d\eta \text{(111)}$$

Cette intégrale double est facile à calculer car les deux variables sont intégrables séparément :

$$\underline{U}(f_x,f_y)=\frac{C}{\lambda z{z}_s}\left[{\int }_{-\frac{a}{2}}^{\frac{a}{2}}exp(-j2\pi f_x\xi )d\xi \right]\left[{\int }_{-\frac{b}{2}}^{\frac{b}{2}}exp(-j2\pi f_y\eta )d\eta \right]\text{(112)}$$

Les deux intégrales sont similaires et on obtient :

$$\underline{U}(f_x,f_y)=\frac{C}{\lambda z{z}_s}ab\frac{sin\left(\pi f_{x}a\right)}{\pi f_{x}a}\frac{sin\left(\pi f_{y}b\right)}{\pi f_{y}b}\text{(113)}$$

L'intensité de la lumière est donc :

$$I(f_x,f_y)=\left(\frac{Cab}{\lambda z{z}_s}{)}^2\right(\frac{sin\left(\pi f_{x}a\right)}{\pi f_{x}a}{)}^2(\frac{sin\left(\pi f_{y}b\right)}{\pi f_{y}b}{)}^2\text{(114)}$$

Ce résultat fait intervenir la fonction sinus cardinal, définie par :

$$sinc\left(u\right)=\frac{sinu}{u}\text{(115)}$$

Voici le tracé de cette fonction. La fonction numpy.sinc(x) calcule en fait le sinus cardinal de $\pi x$.

x = np.linspace(-5,5,1000)
figure(figsize=(8,6))
plot(x,np.sinc(x))
grid()
xlabel(r"$\frac{u}{\pi}$",fontsize=16)
ylabel(r"${\rm sinc}(u)$",fontsize=16)
                 

fig16.pdf

Voici la courbe d'intensité en fonction de f_xa, pour f_y fixé :

x = np.linspace(-5,5,1000)
figure(figsize=(8,6))
plot(x,(np.sinc(x))**2)
grid()
xlabel(r"$f_xa$",fontsize=16)
ylabel(r"$I$",fontsize=16)
                 

fig17.pdf

La majeure partie de l'énergie est concentrée dans la tache centrale de diffraction définie par la condition :

$$-1<f_{x}a<1\text{(116)}$$

soit encore :

$$-\frac{\lambda }{a}<\alpha <\frac{\lambda }{a}\text{(117)}$$

Il est intéressant de reprendre le calcul de l'amplitude diffractée si l'ouverture est décalée dans la direction (Ox) d'une distance d :

$$\underline{U}(f_x,f_y)=\frac{C}{\lambda z{z}_s}{\int }_{d-\frac{a}{2}}^{d+\frac{a}{2}}{\int }_{-\frac{b}{2}}^{\frac{b}{2}}exp(-j2\pi (f_x\xi +f_y\eta ))d\xi d\eta \text{(118)}$$

Avec le changement de variable $\xi =d+\xi \text{'}$; l'intégrale s'écrit :

$$\underline{U}(f_x,f_y)=\frac{C}{\lambda z{z}_s}exp(-j2\pi f_{x}d){\int }_{-\frac{a}{2}}^{\frac{a}{2}}{\int }_{-\frac{b}{2}}^{\frac{b}{2}}exp(-j2\pi (f_x\xi \text{'}+f_y\eta ))d\xi \text{'}d\eta \text{(119)}$$

Le facteur exponentiel qui apparaît en facteur n'a pas d'effet sur l'intensité : dans le cadre de l'approximation de Fraunhofer, la figure de diffraction est insensible à la position de l'ouverture sur le plan (0xy). Remarquons que cette propriété n'est pas valable dans le cas de l'approximation de Fresnel.