[  accueil page perso P. Bénard ] [  accueil Site Officiel du Pic de Saint-Barthélemy ] [ Lever de Soleil au Pic de Saint-Barthélemy) ]

La réfraction atmosphérique en montagne

Pierre Bénard

16 Janvier 2006

Cette page montre comment calculer les effets de la réfraction atmosphérique pour connaître la déformation du paysage en montagne. La méthode décrite ici utilise une discrétisation spatiale du rayon lumineux le long de sa trajectoire, mais elle retient la valeur exacte de l'indice de réfraction en chaque point de l'espace, ce qui conduit à des résultats beaucoup plus précis que des méthodes pour lesquelles le champ d'indice de réfraction est décrit par une approximation discrète fixée a priori. La différence entre la méthode utilisée ici et ces méthodes moins précises est abordée dans l'Appendice B.

1 Introduction

Les effets de la réfraction atmosphérique (mirages, déformations, déviations...) sont dus au fait que l'indice de réfraction de l'atmosphère (noté "$n$" dans la suite) n'est pas constant et uniforme dans l'espace, mais présente des variations locales qui viennent dévier les rayons lumineux par rapport à une trajectoire rectiligne. L'indice de réfraction n'est autre que le rapport de la vitesse de la lumière dans le vide à celle observée dans le milieu réfringent (ici, l'atmosphère). Au cours de leur voyage d'un point à un autre, les rayons lumineux adoptent le trajet qui minimise le temps de parcours, et qui n'est donc pas nécessairement la ligne droite si l'indice de réfraction varie dans l'espace dans la région du trajet.

Il est important de souligner qu'il n'existe pas d'expression simple et universelle donnant la valeur exacte de l'indice de réfraction atmosphérique en un point donné et à un instant donné, car l'indice de réfraction dépend d'une multitude de paramètres eux-mêmes variables dans le temps et dans l'espace (température, pression,...). En conséquence, selon l'application que l'on veut faire des calculs de réfraction, il faut bien analyser quelle sont les dépendances de $n$ par rapport aux différents paramètres qui doivent être retenues et quelles sont celles que l'on peut négliger.

Par exemple, selon que l'on voudra comprendre et calculer de mirages inversés, ou bien l'apparition du "rayon vert", ou encore, comme ici, la déformation du paysage en montagne, les paramètres à prendre en considération dans le calcul de l'indice de réfraction ne seront pas les mêmes. Pour calculer le "rayon vert", il faut évidemment tenir compte des variations de l'indice de réfraction en fonction de la longueur d'onde de la lumière, alors que cet aspect peut être négligé pour calculer la déformation du paysage en montagne, et ainsi de suite...

La formule donnant l'indice de réfraction que nous utilisons comme point de départ pour les calculs ci-dessous ne saurait donc être prise pour argent comptant en vue d'une application différente de celle qui en est faite ici. Néanmoins, en vue d'applications différentes, on pourra toujours tenter de s'inspirer des considérations conduisant à la justification de la formule utilisée ici, et qui sont exposées tout en bas de ce document, dans l'appendice A.

Avec une approximation raisonnable pour ce que nous voulons en faire ici, nous considérerons que les variations de l'indice de réfraction sont uniquement dues à la variation de la densité de l'air sec dans l'atmosphère, et nous supposerons en outre que cette densité est celle de "l'atmosphère standard", c'est à dire une atmosphère homogène horizontalement, et présentant un profil vertical proche de celui celui qui est observé en moyenne sur notre planète.

Ce faisant, nous négligeons donc les effets de réfraction liés à la présence d'autres constituants que l'air sec dans l'atmosphère, notamment la vapeur d'eau. Nous négligeons également de prendre en compte l'état thermique exact de l'atmosphère pour nous concentrer sur l'état standard (qui est un état moyen). Donc l'effet calculé peut différer légèrement d'une observation particulière si l'observation se déroule alors que l'état de l'atmosphère diffère notablement de l'état standard. En outre, comme évoqué plus haut, nous négligeons la dépendance de l'indice de réfraction par rapport à la longueur d'onde du rayonnement.

2 Données du problème

Pour calculer les effets de la réfraction atmosphérique, il faut d'abord connaître comment varie l'indice de réfraction de l'air "$n$" en fonction du lieu, puis il faut connaître la loi de propagation de la lumière dans l'air. Ces deux aspects sont abordés dans cette section.

2.1 Expression de l'indice de réfraction atmosphérique

Pour les problèmes qui nous intéressent, comme expliqué dans l'Appendice A, l'indice de réfraction peut être estimé par la formule suivante:

$$n(z)-1 = (n_{00} -1)\frac{\rho_0}{\rho_{00}} \left(1 - \frac{\gamma z}{T_0} \right)^{(g/ \gamma R)-1} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {\rm si \;} z \leq z_T$$

$$n(z)-1 = (n_{00} -1) \frac{\rho_T}{\rho_{00}}\; \exp \left[- \frac{g }{R T_T}(z-z_T) \right]\;\;\;\;\;\;\;\;\;\;\;\;\; {\rm si \;} z > z_T$$ (1)

où la signification physique des constantes est donnée dans l'Appendice A. Dans cette formule, $z$ est l'altitude (en m), et les températures $T$ sont des températures absolues, en Kelvin.

Pour les applications, nous choisirons les valeurs suivantes:
- $p_{00}=$101325. Pa (soit 1 Atm)
- $T_{00}=$273.16 K (soit $0^\circ$ C)
- $n_{00}=$1.00029
- $z_T=$10000 m
- $\gamma =$ 0.0065 K/m

Ainsi, les valeurs de $n$ au niveau de la mer sont, avec cette formule:

- $T_0 =$ 273.16 K (soit $0^\circ$ C) : $n$=1.00029

- $T_0 =$ 283.16 K (soit $10^\circ$ C) : $n$=1.00028

- $T_0 =$ 293.16 K (soit $20^\circ$ C) : $n$=1.00027

Ensuite, l'indice diminue rapidement avec l'altitude: pour $T_0=283.16$ K (soit $10^\circ$ C), on a:

- $z$ = 0 km: $n$=1.00028

- $z$ = 1 km: $n$=1.00025

- $z$ = 2 km : $n$=1.00023

- $z$ = 3 km : $n$=1.00021

Lorsqu'on applique la forme (1) indiquée ci-dessus pour l'indice de réfraction, on s'aperçoit (après avoir effectué tous les calculs décrits ci-après) que les résultats dépendent en fait assez peu de la température au niveau de la mer $T_0$. On peut donc utiliser partout $T_0$=283.16 K sans commettre une erreur trop importante, et c'est ce que nous faisons désormais.

2.2 Propagation de la lumière dans l'air

La loi de propagation de la lumière dans l'air s'écrit:

$$\frac{d n {\bf u}}{ds} = {\bf grad} (n)$$ (2)

où $u$ est le vecteur unitaire tangent au rayon de lumière et $s$ est l'abcisse curviligne le long du rayon de lumière. Le symbole ${\bf grad}(n)$ désigne le vecteur gradient de l'indice de réfraction.

Dans notre cas, $n$ ne dépend que de $z$ et donc

$${\bf grad} (n) = \frac{d n}{d z} \; {\bf k}$$ (3)

où ${\bf k}$ est le vecteur unitaire dans la direction verticale ascendante. L'expression de la dérivée verticale de l'indice de réfraction en fonction de l'altitude est donnée par l'Appendice A:

$$\frac{dn}{dz} = \frac{1}{T_0}\left(\gamma - \frac{g}{R} \right)\left(n_{00}-1 \ri... ...\right)^{(g/ \gamma R)-2} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {\rm si \;} z \leq z_T$$

$$\frac{dn}{dz} = \frac{- g}{R T_T}(n_{00} -1) \frac{\rho_T}{\rho_{00}}\; \exp \left[- \frac{g }{R T_T} (z-z_T)\right]\;\;\;\;\;\;\;\;\;\;\;\;\; {\rm si \;} z > z_T$$ (4)

D'autre part, si $P$ est un point courant sur le rayon optique et ${\bf OP}$ le vecteur joignant $O$ à $P$, on a, par construction:

$$\frac{d OP}{d s} = {\bf u}$$ (5)

3 Principes de base du calcul d'un rayon optique

On se place dans le plan formé par le centre de la Terre, le point d'observation $P_0$ et le point observé $P_F$. On définit dans ce plan un repère Cartésien direct $(O,x,y)$ ou $O$ est le centre de la Terre et $Oy$ pointe vers $P_0$. Les vecteurs n'ont donc que deux composantes: selon $Ox$ et selon $Oy$. Les variables du problème sont $n$, ${\bf u}$ et le vecteur position des points sur le rayon optique: ${\bf OP}$. Les vecteurs étant décomposés selon $x$ et $y$, on a: ${\bf u} = (u_x, u_y)$ et ${\bf OP} = (OP_x, OP_y)$.

Pour chaque point, connaissant ses coordonnées $(x,y)$, il faut être capable de connaître son altitude $z$. La terre étant ici supposée sphérique de rayon $a$, l'altitude $z$ d'un point de coordonnées $(x,y)$ est donnée par:

$$z = \sqrt{x^2 + y^2}- a$$ (6)

La Figure 1 montre le positionnement de $x$, $y$ et la signification de $z$.

Figure 1: Coordonnées $(x,y)$ et altitude $z$ pour un point $P$ donné. Le rayon terrestre est noté $a$.

Rappelons que l'angle que forme le rayon d'une visée par rapport au plan horizontal, s'appelle le "site" de cette visée (c'est un terme d'artillerie). Pour les visées dirigées vers le bas, le site est négatif.

Le but du problème est de déterminer le rayon lumineux joignant $P_F$ à $P_0$. Ceci nous permettra de résoudre l'ensemble des problèmes que l'on peut se poser:

• $\bullet$ trouver la distance et le site de l'horizon marin vu du sommet du Pic;
• $\bullet$ trouver le haussement apparent du soleil lorsqu'il se couche, vu du sommet du Pic;
• $\bullet$ trouver le site d'un point donné vu du sommet du Pic, et par comparaison avec le site obtenu en négligeant la réfraction, trouver le haussement apparent de ce point.

Comme les variations de $n$ n'ont pas une expression simple, on ne peut pas espérer résoudre simplement le problème par une méthode analytique (c'est à dire à la main, avec un papier et un crayon). Donc, on adopte une méthode numérique: on découpe le rayon lumineux en une succession de petits trajets de longueur $ds$ et pour chacun de ces petits trajets on applique les équations de propagation grâce à une discrétisation de ces équation sous forme de différences finies, en appliquant les valeurs locales des variables nécessaires au calcul. Pour résoudre l'équation il faut alors un ordinateur, et il faut écrire le programme qui résoudra le problème.

Un point important est de bien veiller respecter à le principe de retour inverse de la lumière, sous peine d'imprécision dans les résultats: on s'impose donc de trouver le même résultat selon que l'on vise du point $P_0$ ou du point $P_F$, car si $P_0$ voit $P_F$, alors $P_F$ doit voir $P_0$ par le même rayon lumineux.

4 Calcul discrétisé d'un rayon optique pour une visée donnée

Ici, on commence par s'attaquer à un problème plus simple, qui est de calculer la trajectoire d'un rayon optique pour une visée donnée. On se fixe donc le point de visée $P_0$ et le site de visée $\beta$, et on calcule d'abord la trajectoire suivie par le rayon lumineux, de proche en proche, comme expliqué ci-dessus. Ensuite, puisqu'il faut bien s'arrêter de calculer au bout d'un certain temps, on spécifie une condition d'arrêt de calcul.

4.1 Trajectoire discrétisée d'un rayon optique pour une visée donnée

Rappelons que les équations à discrétiser sont:

$$\frac{d n {\bf u}}{ds} = {\bf grad}(n)$$ (7)

$$\frac{d {\bf OP}}{ds} = {\bf u}$$ (8)

On divise le rayon optique en petits segments de longueur $\delta s$, et donc on suppose que le rayon est décrit par une suite de points $P_{i}$ et de vecteurs unitaires tangents ${\bf u}_i$, qu'il nous faut calculer. L'indice de réfraction à ce point $P_i$ est noté $n_i$. Les équations discrétisées s'écrivent alors:

$$\frac{n_{i+1} {\bf u}_{i+1} - n_{i} {\bf u}_{i}}{\delta s} = \frac{{\bf grad}(n_{i+1}) +{\bf grad}(n_{i}) }{2}$$ (9)

$$\frac{{\bf OP}_{i+1} - {\bf OP}_{i}}{\delta s} = \frac{{\bf u}_{i+1} +{\bf u}_{i}}{2}$$ (10)

On suppose connus $P_i$ et ${\bf u}_i$ et on veut calculer le point suivant: $P_{i+1}$, et le vecteur unitaire tangent au rayon au point suivant: ${\bf u}_{i+1}$, ainsi que l'indice de réfraction au point suivant $n_{i+1}$. Comme on s'est imposé une discrétisation "implicite" (c'est à dire symétrique) afin d'assurer la contrainte du retour inverse de la lumière, les inconnues apparaissent à plusieurs endroits, et il n'est pas simple de résoudre le système précédent par une méthode directe. On utilise donc une méthode itérative. D'abord on calcule un estimateur (noté $P_{i+1}^{(0)}$) du point suivant $P_{i+1}$ et les variables associées:

• on calcule $P_{i+1}^{(0)} = P_i + {\bf u}_i \; \delta s$;
• on calcule $n_{i+1}^{(0)}$ et ${\bf grad} \left[n_{i+1}^{(0)} \right]$ en $P_{i+1}^{(0)}$;
• on calcule ${\bf u}_{i+1}^{(0)}$ par: $n_{i+1}^{(0)}{\bf u}_{i+1}^{(0)} =n_{i} {\bf u}_{i} + 0.5 \; \delta s \; \left... ...{\bf grad} \left[n_{i+1}^{(0)} \right]+{\bf grad} \left[n_{i} \right]\right \}$;

Une fois connu l'estimateur $P_{i+1}^{(0)}$ et toutes les variables associés, on applique le processus itératif en lui-même. Pour un indice d'itération noté $(k)$, on applique l'algorithme itératif suivant:

• on calcule $P_{i+1}^{(k)} = P_i + 0.5 ( {\bf u}_{i+1}^{(k-1)} + {\bf u}_i ) \; \delta s$;
• on calcule $n_{i+1}^{(k)}$ et ${\bf grad} \left[n_{i+1}^{(k)} \right]$ en $P_{i+1}^{(k)}$;
• on calcule ${\bf u}_{i+1}^{(k)}$ par: $n_{i+1}^{(k)}{\bf u}_{i+1}^{(k)}=n_{i} {\bf u}_{i} + 0.5 \; \delta s \; \left \{ {\bf grad} \left[n_{i+1}^{(k)} \right]+{\bf grad} \left[n_{i} \right]\right \}$;

Du fait de la faible non-linéarité des équations, l'algorithme converge extrêmement rapidement, et en pratique une seule itération est nécessaire pour arriver à une grande précision. Donc, connaissant un point $P_i$ sur le rayon lumineux, le point suivant sur la trajectoire est donné par:

$$P_{i+1} = P_{i+1}^{(1)}$$ (11)

On voit que la contrainte de respect du principe de retour inverse de la lumière n'est qu'imparfaitement respectée puisqu'on ne fait qu'une seule itération, ce qui ne garantit pas que les équations de la discrétisation symétrique sont exactement résolues. Toutefois puisque la convergence du processus itératif est très rapide, la contrainte est respectée avec une bonne précision, en tout cas suffisante pour notre problème.

On voit également que la pécision du calcul dépendra de la longueur des segments $\delta s$ que l'on a choisi. Si on prend $\delta s$ aussi long que la distance entre $P_0$ et $P_F$, la précision sera très mauvaise, puisque le rayon modélisé joignant les deux points ne sera pas une courbe mais un segment de droite. Dans tous les calculs que je décris, j'ai utilisé la valeur $\delta s = 100 m$ ou $200 m$, qui donne une précision largement compatible avec nos objectifs. Bien évidemment, plus on choisit $\delta s$ petit, plus la précision du calcul est bonne, mais aussi, plus le calcul est long, car il faut décomposer le rayon lumineux joignant $P_0$ et $P_F$ en un plus grand nombre de segments, qu'il faut alors calculer un par un à la suite l'un de l'autre.

4.2 Conditions de départ

Pour pouvoir débuter le calcul, il faut spécifier le point de visée qui, par construction, a pour coordonnées:

$$P_0 = (0, a+ z_0)$$ (12)

où $z_0$ est l'altitude du point. Il faut également spécifier l'angle de visée (le site de visée) que nous appelerons $\beta_0$. Les composantes du vecteur ${\bf u}_0$ tangent au rayon au point $P_0$ sont alors:

$${\bf u}_{0x} = \cos \beta_0$$ (13)

$${\bf u}_{0y} = \sin \beta_0$$ (14)

4.3 Condition d'arrêt

Il faut nécessairement spécifier une fin au calcul de la trajectoire du rayon. En pratique, ce qui nous intéresse c'est d'arrêter le calcul une fois que le rayon a franchi la distance qui sépare le point de visée $P_0$ du point observé $P_F$. Seulement comme le rayon est courbé, ce n'est pas la distance le long du rayon qu'il faut considérer, mais la distance absolue sur la Terre. On commence donc par mesurer la distance sur la Terre du point observé $P_F$ par rapport au point de visée $P_0$. En fait c'est la distance angulaire $\alpha_{\rm OF}$ entre ces deux points que l'on utilise, c'est à dire l'angle formé par les deux rayons $(OP_0, OP_F)$, où $O$ est, rappelons le, le centre de la Terre. Si $d_{\rm OF}$ est la distance sur la Terre en les deux points, cette distance angulaire $\alpha_{\rm OF}$ est donnée par:

$$\alpha_{\rm OF} = \frac{d_{\rm OF}}{a}$$ (15)

où $a$ et $d$ sont en km et $\alpha_{\rm OF}$ est en Radians. Si on ne connaît pas la distance terrestre (qui est difficile à mesurer sur une carte pour des grandes distances, à cause des déformations liées à la projection utilisée pour obtenir une carte plate), on peut alors faire un calcul plus compliqué à partir des coordonnées géographiques des points $P_0$ et $P_F$. Mais jusqu'à 500 km, on peut normalement se contenter de mesurer directement la distance avec une règle sur la carte sans trop d'erreur.

Pendant le calcul de la trajectoire de notre rayon, nous pouvons donc calculer pour chaque point $P_i = (x_i,y_i)$ la distance angulaire sur la Terre déjà parcourue par le rayon:

$$\alpha_{i} = {\rm ArcTan} (\frac{x_i}{y_i})$$ (16)

Il suffit donc d'arrêter le calcul dès que $\alpha_i > \alpha_{\rm OF}$. C'est à dire dès que le rayon lumineux aura dépassé la verticale du point observé.

Il est important de remarquer que comme ici nous avons spécifié $\beta_0$ de manière libre ou au hasard, rien ne garantit que notre rayon passe bien au point $P_F$ lorsque la distance $\alpha_{\rm OF}$ est franchie: il peut très bien passer au-dessus si $\beta_0$ a été choisi trop grand, soit au-dessous dans le cas inverse. Dans la prochaine section on va s'intéresser à trouver la valeur de $\beta_0$ qui permet au rayon d'atteindre exactement le point $P_F$.

On peut calculer l'erreur par rapport à la cible $P_F$. Pour cela on enregistre l'altitude $z_{\rm prec}$ du rayon au dernier point de trajectoire précédant la cible $P_F$, et dont la distance angulaire avec $P_0$ est notée $\alpha_{\rm prec}$. De même on enregistre l'altitude du rayon au point suivant $z_{\rm suiv}$, qui est le premier point de trajectoire ayant dépassé la cible, et dont la distance angulaire avec $P_0$ est notée $\alpha_{\rm suiv}$. On calcule alors l'altitude du rayon sur ce segment au moment où il dépasse la cible par une interpolation linéaire:

$$z= z_{\rm prec} + \frac{ (\alpha_{\rm OF} - \alpha_{\rm pre... ...{\rm suiv} - \alpha_{\rm prec})}(z_{\rm suiv}-z_{ \rm prec})$$ (17)

5 Calcul du rayon joignant les deux points désirés

Maintenant on calcule la trajectoire du rayon qui joint exactement les deux points qui nous intéressent: $P_0$ et $P_F$. Pour cela nous adoptons simplement une méthode par "essai-erreur", en réduisant successivement l'erreur au moyen d'un nouveau processus itératif. Rappelons que l'on note $z_0$ l'altitude du point d'observation $P_0$, et $z_1$ l'altitude du point cible $P_F$.

Tout d'abord, on spécifie un $\beta_{\rm max}$ qui vise à coup sûr au-dessus de la cible, et un $\beta_{\rm min}$ qui vise à coup sûr au-dessous de la cible. En pratique, on prend par exemple $\beta_{\rm max}=0.05$ Rad, et $\beta_{\rm min}= -0.05$ Rad. Pour ces deux valeurs initiales, on calcule l'altitude obtenue lorsque le rayon atteint la distance de la cible, que l'on note $z_{\rm min}$ et $z_{\rm max}$. On en déduit donc les erreurs survenues pour ces deux visées.

Ensuite, on applique le processus itératif d'essais-erreurs lui-même:

• $\bullet$ on calcule le meilleur site possible pour notre cible à partir des informations connues, c'est à dire: $\beta= \beta_{\rm max} - (\beta_{\rm max}-\beta_{\rm min})(z_{\rm max}-z_1)/(z_{\rm max}-z_{\rm min})$;
• $\bullet$ on calcule l'altitude $z$ obtenue à la distance angulaire $\alpha_{\rm OF}$ avec ce site;
• $\bullet$ si $z>z1$ on affecte à $\beta_{\rm max}$ la nouvelle valeur $\beta$, sinon on affecte a $\beta_{\rm min}$ la nouvelle valeur $\beta$.

On réitère ce processus un nombre fixé de fois (en pratique une dizaine d'itérations suffisent). En imprimant l'erreur d'altitude finale par rapport à la cible, on s'assure que la précision atteinte est celle que l'on désire, sinon, il faut augmenter le nombre d'itérations. L'erreur diminue très vite, d'un facteur 10 à 100 à chaque nouvelle itération. On notera que cette méthode est exactement celle qui est utilisée par les artilleurs, dont le problème se pose d'ailleurs en termes très proches (cette méthode nécessite toutefois la possibilité de pouvoir estimer visuellement l'erreur du tir d'artillerie par rapport à la cible, au moyen d'observateurs avancés, ou placés sur des hauteurs avoisinantes).

Nous sommes maintenant capables de calculer le site de visée $\beta$ qu'il faut adopter à partir de notre point d'observation pour que le rayon lumineux atteigne effectivement le point observé.

6 Applications

6.1 Site apparent d'un point donné sur la Terre.

Maintenant, nous savons déterminer la trajectoire d'un rayon lumineux pour qu'il aboutisse en un point donné, et donc nous pouvons calculer le site de visée $\beta$ pour observer ce point: il suffit de spécifier l'altitude $z_0$ du point d'observation (pour le Pic de Saint-Barthélemy, 2348m), la distance angulaire du point observé $\alpha_D$ et son altitude $z_1$. L'algorithme que nous avons détaillé ci-dessus donne alors le site apparent du point observé.

6.2 Distance de l'horizon pour une altitude donnée.

Si on veut calculer la distance de l'horizon optique observé à partir d'une altitude donnée $z_{00}$, on peut se servir du même programme en le modifiant très légèrement: On suppose une altitude du point de visée $z_0$ égale à zéro, et un site de visée $\beta$ égal à zéro lui aussi. On peut calculer la trajectoire du rayon de visée par la méthode indiquée précédemment. Puisque l'observateur est au ras de la mer et vise à l'horizontale, le rayon de visée va s'élever peu à peu par rapport au sol. Il suffit alors de détecter à quelle distance le rayon passe exactement à l'altitude souhaitée $z_{00}$. Alors, en vertu du principe de retour inverse de la lumière, cette distance est aussi la distance apparente de l'horizon vu du point d'altitude $z_{00}$. Ensuite, reprenant le programme initial, en fixant pour le point d'observation l'altitude cible $z_0 = z_{00}$ et pour le point cible la distance obtenue précédemment et l'altitude $z_1=0$, on peut calculer le site apparent de l'horizon optique, vu du point d'observation à l'altitude $z_{00}$.

Si l'observateur est très près du niveau de la mer (quelques centimètres par exemple) l'horizon est situé dans une direction de visée pratiquement horizontale et extrêmement proche de lui. Dans ce cas, les rayons lumineux considérés sont tellement courts que la réfraction joue un rôle pratiquement négligeable. En revanche, pour des altitudes supérieures à 1000m, la différence avec la distance et le site que l'on trouverait de manière purement géométrique (c.à.d. sans tenir compte de la réfraction) commence à être assez significative. Ainsi, pour le Pic de Saint-Barthélemy, l'horizon géométrique est à 173 km avec un site de $-1.552$ degrés, tandis que l'horizon optique se trouve repoussé à 190.3 km, avec un site $\beta_h = -1.422$ degrés. La distance angulaire de l'horizon optique est alors $\alpha_D=1.703$ degrés. On voit que l'horizon marin vu du sommet du Pic de Saint-Barthélemy est rehaussé de $0.13$ degrés environ à cause de la réfraction atmosphérique.

6.3 Haussement du soleil à son lever en plaine.

Si dans la méthode de détermination de l'horizon que nous venons juste d'aborder, on fixe une altitude $z_{00}$ supérieure à 20 km environ pour la cible, on s'aperçoit que lorsqu'il a atteint une certaine altitude, le rayon lumineux, devient pratiquement rectiligne, c'est à dire que le vecteur ${\bf u}$ devient pratiquement constant. Ceci est dû au fait qu'à partir d'une certaine altitude, il n'y a plus assez d'air pour dévier les rayons, et c'est alors la loi de propagation (rectiligne) dans le vide qui s'applique.

En calculant l'angle que fait le vecteur ${\bf u}$, lorsqu'il devient constant, avec l'axe $Ox$ (qui est, rappelons-le, la direction de visée dans le problème précédent), on en déduit le haussement apparent d'un objet extra-atmosphérique situé exactement à l'horizon pour l'observateur. C'est donc aussi l'angle de haussement du soleil lorsqu'il se couche ou se lève en mer ou en plaine.

La Figure 2 illustre la situation examinée: l'observateur se trouve au point $H$ (en plaine ou au ras de la mer), et observe l'horizon, donc à l'horizontale, dans la direction $HC$. Toutefois, à cause de la réfraction atmosphérique, le rayons lumineux qu'il voit n'est pas la droite $HC$, mais la courbe en pointillé qui passe par $H$ et $S$. A partir d'une certaine altitude, cette courbe devient droite comme on le voit sur la figure: à partir du point $S$, le rayon suit la direction de la droite $DS$. L'angle (noté $h_s$ sur la figure) de la déviation du rayon entre sa direction initiale $HS$ et sa direction finale $DS$, est le haussement du soleil à son lever en plaine ou en mer.

Figure 2: Différentes notations d'angles pour une observation en mer (au point $H$) ou au sommet du Pic (au point $P$).

Le calcul donne un haussement apparent de $h_s = 0.562$ degrés. C'est à dire que lorsque l'on observe le début du coucher de soleil, celui-ci est en fait déjà "couché" du point de vue géométrique, puisque son diamètre apparent est d'environ 0.5 degrés. Ce résultat est un résultat classique de l'application de la théorie de la réfraction atmosphérique.

6.4 Haussement du soleil à son lever au sommet du Pic.

Pour trouver le haussement du soleil au moment de son lever vu du sommet du Pic de Saint-Barthélemy, il suffit d'ajouter le haussement du soleil au lever en plaine, $h_s$, obtenu ci-dessus , et la déviation du rayon lumineux entre le sommet du Pic et l'horizon optique vu du sommet du Pic.

La situation examinée est encore illustrée sur la Figure 2. Cette fois, l'observateur est en haut du Pic, au point $P$. Il regarde dans la direction de l'horizon optique, c'est à dire dans la direction $PB$ sur la figure. Le rayon lumineux, en pointillé, passe alors par l'horizon optique, qui est le point $H$, puis continue son trajet vers le point $S$, où il devient une droite. Le haussement du soleil à son lever au Pic est donc l'angle entre la droite $PB$ et la droite $DS$. Pour connaître cet angle, il faut connaître l'angle entre $PB$ et la droite $HC$. On voit sur la figure que cet angle entre $PB$ et $HC$ est égal à $\alpha_D -\left \vert \beta_h \right \vert$, c'est à dire la distance angulaire terrestre entre le Pic et l'horizon optique moins la valeur absolue du site de l'horizon optique vu du Pic (notons que cette formule n'est valable que si $\beta_h$ est négatif, c'est à dire si l'horizon est situé à une altitude inférieure au point de visée). La droite $PA$ sur la figure est la droite horizontale passant au point $P$. L'angle $(OP, PA)$ est un angle droit, tout comme l'angle $(OH, HC)$, ce qui explique que l'angle $(PA, HC)$ vaut bien $\alpha_D$.

Le haussement du soleil à son lever vu du Pic est donc égal à:

$$h_P = h_s +\alpha_D - \left \vert \beta_h \right \vert$$ (18)

Le calcul donne un haussement apparent de $h_P = 0.843$ degrés. Le fait que l'observateur soit en altitude augmente encore le haussement du soleil à son lever par rapport à une observation effectuée en plaine. En conséquence, lorsque l'on voit le coucher de soleil vu du sommet du Pic, le soleil est en fait "couché depuis un certain temps" (du point de vue géométrique).

7 Appendice A: Expression de l'indice de réfraction

Comme indiqué dans l'introduction, pour ce que nous voulons en faire ici, nous considérerons que les variations de l'indice de réfraction sont uniquement dues à la variation de la densité de l'air sec dans l'atmosphère, et nous supposerons en outre que cette densité est celle de "l'atmosphère standard", c'est à dire une atmosphère homogène horizontalement, et présentant un profil vertical proche de celui celui qui est observé en moyenne sur notre planète. Il faut donc définir précisément le profil standard utilisé, puis en déduire la formule de variation de l'indice de réfraction en fonction du lieu.

7.1 Indice de réfraction de l'air sec

Dans le vide, l'indice de réfraction vaut 1 par définition. Dans l'air, la vitesse est ralentie et l'indice $n$ est donc supérieur à 1. Des mesures précises on montré que pour une pression $p_{00}$ de 101325. Pa, et une température $T_{00}$ de 273.16 K, l'indice de réfraction de l'air sec $n_{00}$ est donné par:

$$n(\sigma) = 1+ 10^{-8} \left(6432.8 + \frac{2949810}{146-\sigma^2} + \frac{25540}{41-\sigma^2 } \right)$$ (19)

où $\sigma$ est l'inverse de la longueur d'onde exprimée en $\mu$m, c'est à dire:

$$\sigma = [\lambda(\mu{\rm m}) ]^{-1}$$ (20)

Pour le problème qui nous intéresse, on prendra une longueur d'onde $\lambda$ de 0.5 $\mu$m, soit $\sigma=2$. Pour les conditions $(p_{00},T_{00})$, l'indice de réfraction est alors:

$$n_{00} = 1.00029$$ (21)

Ensuite, la propagation de la lumière dans l'air est telle que le ralentissement de la lumière par rapport à la propagation dans le vide est proportionnel à la densité de l'air traversé. Si on note $\rho$ la densité de l'air et $\rho_{00}$ sa densité pour les conditions $(p_{00},T_{00})$, on a:

$$n(\rho)-1 = (n_{00}-1) \frac{\rho}{\rho_{00}}$$ (22)

Enfin, comme l'air est considéré ici comme un gaz parfait, on a:

$$p = \rho R T$$ (23)

où $R$ est la constante des gaz parfaits pour l'air sec ($R \approx$ 287. J/K). On en déduit donc finalement:

$$n(p,T)-1 = (n_{00}-1) \frac{p T_{00}}{T p_{00}}$$ (24)

7.2 Atmosphère standard

Il est important de souligner que nous souhaitons non seulement calculer la déformation du paysage, mais aussi la trajectoire de rayons issus du soleil et parvenant au sol, afin de calculer leur déviation apparente. Cela nécessite donc de connaître l'indice de réfraction pour des couches atmosphériques situées non seulement près du sol, mais aussi très haut au dessus du sol (typiquement jusqu'à un cinquantaine de kilomètres). Pour cela, nous sommes alors obligés de connaître l'indice de réfraction et donc la densité de l'air jusqu'à ces altitudes élevées.

Nous notons $z$ l'altitude au-dessus du niveau de la mer.

Nous définirons notre atmosphère standard par une atmosphère d'air sec (considéré comme un gaz parfait) en équilibre hydrostatique et composée de deux couches superposées et de caractéristiques thermiques différentes, dont l'interface (appelée tropopause) est située à l'altitude $z_T$. On a alors

• du sol jusqu'à l'altitude $z_T$, la première couche, supposée représenter la troposphère, et dans laquelle le gradient vertical de température sera pris constant et vaudra $\gamma =$ -6.5 K/km. Donc, quand on s'élève de 1 km, la température y diminue de 6.5 degrés. Cela correspond à la moyenne planétaire de la variation thermique avec l'altitude.
• au-dessus de l'altitude $z_T$, la seconde couche supposée représenter la stratosphère, dans laquelle le gradient vertical de température sera considéré comme nul. La stratosphère est donc ici considérée comme une couche isotherme, ce qui correspond aussi grosso-modo aux observations moyennes.

Les subtilités liées aux couches supérieures de l'atmosphère (mésosphère, où la température se met à croître avec l'altitude, puis ionosphère, où l'atmosphère devient ionisée et où l'indice de réfraction se met à dépendre surtout de la densité d'électrons) sont ignorées, car elles affectent peu les phénomènes que nous étudions. On considère donc qu'à partir de l'altitude $z_T$, l'atmosphère est isotherme jusqu'à l'infini. Notons $T_T$ la température au niveau $z=z_T$, et $p_T$ la pression au niveau $z=z_T$.

Remarque: Si l'on n'est pas intéressé par les calculs de position apparente du soleil, et que l'on s'intéresse seulement aux calculs de déformation du paysage, alors on peut ignorer la couche isotherme (stratosphère) dans les calculs. En effet, les reliefs terrestres sont trop peu élevés pour pénétrer jusqu'à la stratosphère, et donc les rayons lumineux d'un point à un autre du paysage sont toujours confinés dans la troposphère.

Nous appellerons $p_0$ la pression au niveau de la mer et $T_0$ la température au niveau de la mer. La valeur de ces deux quantités n'est pas spécifiée pour l'instant afin de pouvoir les choisir proches des conditions observées.

Notre atmosphère standard est donc définie par le profil thermique $T(z)$ suivant:

$$T(z) = T_0 - \gamma z \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {\rm si \;} z \leq z_T$$

$$T(z) = T_T = T_0 - \gamma z_T \;\;\;\;\; {\rm si \;} z > z_T$$ (25)

Il faut maintenant calculer le profil barique de l'atmopshère. Pour cela, nous utilisons le fait que nous avons supposé l'atmopshère en équilibre hydrostatique (par exemple au repos). On a alors:

$$\frac{d p}{d z} = - \rho(z) g$$ (26)

où $g$ est l'accélération de la gravité $g=9.81$ m/s/s.

Dans la troposphère, (26) donne:

$$\frac{d \ln p}{d z} = - \frac{g}{R (T_0 - \gamma z)}$$ (27)

et donc, puisque $\gamma$ est non nul:

$$p(z) = p_0 \left(1 - \frac{\gamma z}{T_0} \right)^{(g/ \gamma R)}$$ (28)

Dans la stratosphère isotherme, (26) donne:

$$\frac{d \ln p}{d z} = - \frac{g}{R T_T}$$ (29)

et donc:

$$p(z) = p_T \; \exp \left[- \frac{g }{R T_T} (z-z_T)\right]$$ (30)

où $p_T = p_0 \left(1 - \frac{\gamma z_T}{T_0} \right)^{(g/ \gamma R)}$ est déterminé en substituant $z=z_T$ dans (28), car le profil de pression doit être continu en $z=z_T$. Le profil de pression est donc donné par:

$$p(z) = p_0 \left(1 - \frac{\gamma z}{T_0} \right)^{(g/ \gamma R)} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {\rm si \;} z \leq z_T$$

$$p(z) = p_T \; \exp \left[- \frac{g }{R T_T} (z-z_T)\right]\;\;\;\;\;\;\;\;\; {\rm si \;} z > z_T$$ (31)

Finalement, le profil de densité s'exprime alors par:

$$\rho(z) = \rho_0 \left(1 - \frac{\gamma z}{T_0} \right)^{(g/ \gamma R)-1} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {\rm si \;} z \leq z_T$$

$$\rho(z) = \rho_T \; \exp \left[- \frac{g }{R T_T}(z-z_T) \right]\;\;\;\;\;\;\;\;\;\;\;\;\; {\rm si \;} z > z_T$$ (32)

où $\rho_T = \rho_0 \left(1 - \frac{\gamma z_T}{T_0} \right)^{(g/ \gamma R)-1}$ est la densité au niveau de la tropopause.

7.3 Indice de réfraction de l'air

En utilisant la formule (22), on obtient le profil de l'indice de réfraction atmosphérique:

$$n(z)-1 = (n_{00} -1)\frac{\rho(z)}{\rho_{00}},$$ (33)

soit:

$$n(z)-1 = (n_{00} -1)\frac{\rho_0}{\rho_{00}} \left(1 - \frac{\gamma z}{T_0} \right)^{(g/ \gamma R)-1} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {\rm si \;} z \leq z_T$$

$$n(z)-1 = (n_{00} -1) \frac{\rho_T}{\rho_{00}}\; \exp \left[- \frac{g }{R T_T}(z-z_T) \right]\;\;\;\;\;\;\;\;\;\;\;\;\; {\rm si \;} z > z_T$$ (34)

Le gradient vertical de l'indice de réfraction s'écrit alors:

$$\frac{dn}{dz} = \frac{1}{T_0}\left(\gamma - \frac{g}{R} \right)\left(n_{00}-1 \ri... ...\right)^{(g/ \gamma R)-2} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {\rm si \;} z \leq z_T$$

$$\frac{dn}{dz} = \frac{- g}{R T_T}(n_{00} -1) \frac{\rho_T}{\rho_{00}}\; \exp \left[- \frac{g }{R T_T} (z-z_T)\right]\;\;\;\;\;\;\;\;\;\;\;\;\; {\rm si \;} z > z_T$$ (35)

8 Appendice B: Comparaison avec d'autres méthodes de calcul

La méthode de calcul utilisée ici n'est pas la seule possible: il existe d'autres méthodes pour évaluer la trajectoire de rayons lumineux dans des milieux réfringents. Il existe notamment une méthode légèrement plus simple que celle décrite ci-dessus, et qui consiste à considérer l'atmosphère comme une série de $N$ couches superposées, à l'intérieur desquelles l'indice est uniforme. Nous appellerons cette méthode "méthode des dioptres superposés", car à l'interface de chaque couche d'indice uniforme, il existe une discontinuité de la valeur de l'indice, ce qui crée donc un dioptre. Dans cette méthode, l'atmosphère est donc considérée comme une superposition de $N$ dioptres d'indice uniforme..

La méthode de calcul que nous avons décrit illustrée ici est beaucoup plus précise et plus universelle que celle des dioptres superposés, car elle peut s'appliquer telle quelle dès que l'expression exacte de l'indice de réfraction est connue. Cette méthode est notamment valable quelle que soit la distance entre l'observateur et la cible observée et quelle que soit la direction du rayon calculé.

La méthode des dioptres superposés est moins universelle en ce sens par exemple qu'elle ne permet pas de décrire correctement la trajectoire d'un rayon lumineux qui demeure dans une même couche tout au long de son parcours. Par construction, avec cette méthode, la trajectoire d'un tel rayon sera un segment de droite, alors que dans la réalité, du fait des variations spatiales continues de l'indice de réfraction, la trajectoire du rayon est incurvée. Cette méthode des dioptres superposés est donc plutôt mieux adaptée au calcul des rayons inclinés que pour le calcul des rayons proches de l'horizontale.

Or, dans la réalité, les rayons qui ont la plus forte courbure locale sont ceux qui sont les plus proches de l'horizontale. Ceci est démontré sur le lien suivant:

  courbure locale d'un rayon lumineux dans un milieu à réfringence continument variable 
La méthode des dioptres superposés ne permet pas de rendre compte de cette réalité.

On peut dire que dans la méthode des dioptres superposés, la discrétisation, étant adoptée a priori et selon la verticale, n'est pas pertinente pour des rayons qui ont des trajets horizontaux, car ces rayons sont alors très mal échantillonnés. Par exemple si on adopte des couches de 100m d'épaisseur, des rayons horizontaux proches de la surface terrestre peuvent voyager, avec ce modèle, sur une distance de plus de 140 km sans être aucunement échantillonnés: l'échantillonage est alors très médiocre pour ces rayons, à l'endroit même où il devrait être le meilleur possible, du fait de la courbure locale maximale de ces rayons dans la réalité, comme on vient de le voir.

La méthode décrite dans la note ci-dessus assure en revanche un échantillonage uniforme des rayons sur toute leur longueur, ce qui conduit à une méthode de discrétisation de problème bien mieux posée et bien plus pertinente.

Ajoutons que d'un point de vue conceptuel, cette méthode des dioptres superposés pourrait laisser croire que les rayons qui se propagent dans une direction proche de l'horizontale sont moins courbés que ceux qui ont une inclinaison notable, alors que c'est précisément l'inverse qui est vrai. En fait, cette méthode étant imprécise par nature, c'est son imprécision qui produit ce résultat. Alors, si on se contente d'examiner les résultats produits par cette méthode et qu'on les adopte comme vérité, on en tirera la conclusion que les rayons horizontaux sont moins courbés que leurs homologues inclinés, ce qui est faux en réalité.

En revanche, pour des trajets de rayons dans l'atmosphère suffisament longs et inclinés par rapport à l'horizontale, cette méthode des dioptres superposés, bien qu'imprécise, donne des résultats relativement satisfaisants en première approximation, car alors la discrétisation est pertinente sur une partie importante du trajet, du fait de l'inclinaison des rayons. Pour de tels rayons, la majeure partie du trajet s'effectue de façon inclinée et donc le modèle des dioptres superposés n'est pas trop imprécis, pourvu que le nombre de couches $N$ soit choisi de façon adéquate par rapport au problème.

Reste qu'avec le développement des moyens automatiques de calculs, il n'y a plus de raison de continuer à utiliser la méthode des dioptres superposés, potentiellement mal échantillonée et imprécise, comme on vient de le voir. Cette méthode était bien adaptée pour des calculs "à la main" ou avec les très faibles moyens de calculs dont on disposait avant l'ère de l'électronique, du fait du nombre très modique d'opérations qu'elle requiert, mais elle nécessitait de grandes précautions dans son application. De nos jours, des méthodes légèrement plus couteuses en nombre d'opérations mais plus universelles et plus précises devraient être logiquement privilégiées.