Dioptres centrés : optique géométrique

2.a. Translation et réfraction

On considère un dioptre sphérique de courbure C (nulle pour un dioptre plan), séparant deux milieux d'indice de réfraction n₁ et n₂. On note S le sommet, c'est-à-dire l'intersection du dioptre avec l'axe optique Oz, et C son centre de courbure.

Figure pleine page

Par convention, l'indice n₁ désigne l'indice du milieu d'où vient le rayon, éventuellement à droite du dioptre.

La courbure est définie de manière algébrique, positive lorsque le centre est à droite du sommet, cela indépendamment du sens de la lumière (qui peut venir de la droite) :

$$C=\frac{1}{R}=\frac{1}{\overline{SC}}$$

Soit un rayon lumineux partant d'un point P₁(x₁,y₁,z₁) dans une direction donnée par le vecteur unitaire U₁(u₁,v₁,w₁). Ce rayon rencontre la surface du dioptre au point P₂(x₂,y₂,z₂) où la normale unitaire est N, avec un angle d'incidence θ₁. Le rayon réfracté est défini par le vecteur unitaire U₂(u₂,v₂,w₂).

On recherche tout d'abord les coordonnées de P₂ sous la forme :

$$\begin{array}{cc}\hfill & x_2=x_1+t{u}_1\hfill \\ \hfill & y_2=y_1+t{v}_1\hfill \\ \hfill & z_2=z_1+t{w}_1\hfill \end{array}$$

où t est la distance parcourue de P₁ à P₂. Cette distance est obtenue en considérant l'intersection du rayon avec la surface du dioptre.

Considérons tout d'abord le cas d'un dioptre plan, pour lequel la courbure est nulle. L'equation du dioptre est z₂=z_s et la distance est donc:

$$t=\frac{z_s-z_1}{w1}$$

La valeur doit être positive, sinon le rayon de coupe pas le dioptre.

Lorsque la courbure est non nulle, il faut considérer l'intersection du rayon avec la surface sphérique d'équation :

$$x_2^2+y_2^2+(z_2-z_S-R{)}^2=R^2$$

où z_S définit la position du sommet. On obtient ainsi une équation du second degré vérifiée par t:

$$t^2+2bt+d=0$$

Les coefficients de cette équation sont :

$$\begin{array}{cc}\hfill & b=x_1{u}_1+y_1{v}_1+(z_1-z_S)w_1-R{w}_1\hfill \\ \hfill & d=x_1^2+y_1^2+{(z_1-z_S)}^2-2R(z_1-z_S)\hfill \end{array}$$

Le discriminant de l'équation est :

$$\Delta =b^2-d$$

On se place dans le cas où Δ est positif, sinon le rayon incident ne rencontre pas le dioptre. Les deux solutions sont :

$$t=-b\pm \sqrt{b^2-d}$$

On exclue le cas d'une racine double, qui correspond à un rayon tangent à la sphère. Le choix entre les deux solutions se fait sachant que t doit être positif, et que l'intersection doit se situer sur l'hémisphère où se trouve le sommet. Cette dernière condition s'écrit :

$$(z_2-z_c)(z_S-z_c)>0$$

Si deux solutions positives vérifient cette condition, il faut choisir la plus petite.

Voyons à présent le calcul du vecteur unitaire U₂ définissant le rayon réfracté, obtenu en appliquant les lois de Descartes. On définit tout d'abord le vecteur normal unitaire par :

$$\overrightarrow{N}=\frac{\overrightarrow{C{P}_2}}{R}$$

Ce vecteur est orienté dans le sens de z décroissant, quelque soit le signe de la courbure. Le rayon réfracté est dans le plan d'incidence, donc il existe un nombre réel K tel que :

$$K\overrightarrow{N}=n_1\overrightarrow{U_1}-n_2\overrightarrow{U_2}$$

Par projection sur la normale, on obtient K en fonction des angles d'incidence et de réfraction :

$$K=-n_{1}cos{\theta }_1+n_{2}cos{\theta }_2$$

Le cosinus de l'angle d'incidence est calculé avec le produit scalaire suivant :

$$-\overrightarrow{U_1}\cdot \overrightarrow{N}=cos{\theta }_1=w_1\frac{R+z_S-z_2}{R}-u_1\frac{x_2}{R}-v_1\frac{y_2}{R}$$

Les deux expressions précédentes sont valables si le rayon incident se propage vers la droite (w₁>0). Dans le cas contraire, les expressions opposées doivent être utilisées.

Le cosinus de l'angle de réfraction est obtenu par

$$\begin{array}{cc}\hfill & sin{\theta }_2=\frac{n_1}{n_2}sin{\theta }_1\hfill \\ \hfill & cos{\theta }_2=\sqrt{1-{sin}^2{\theta }_2}\hfill \end{array}$$

puis on en déduit K avec la relation obtenue plus haut. Finalement les composantes du vecteur unitaire du rayon réfracté sont données pour un dioptre sphérique par :

$$\begin{array}{cc}\hfill n_2{u}_2& =n_1{u}_1-CK{x}_2\hfill \\ \hfill n_2{v}_2& =n_1{v}_1-CK{y}_2\hfill \\ \hfill n_2{w}_2& =n_1{w}_1-CK(z_2-R-z_S)\hfill \\ \hfill & =n_1{w}_1-CK(z_2-z_s)+K\hfill \end{array}$$

Pour finir, on remarque que la réflexion par un miroir sphérique se traite avec les mêmes équations en posant :

$$n_2=-n_1$$

Par exemple, si la lumière arrive dans le sens positif (z croissant) sur un miroir placé dans l'air, on posera :

$$\begin{array}{cc}\hfill & n_1=1\hfill \\ \hfill & n_2=-1\hfill \end{array}$$

Si la lumière arrive sur le miroir dans le sens négatif, les indices sont opposés. L'indice négatif est toujours associé à la lumière se propageant dans le sens négatif. Les rayons de courbure des dioptres étant définis indépendamment du sens de la lumière, la courbure d'un dioptre donné change de signe avec le sens de la lumière.

2.b. Approximation paraxiale

L'approximation paraxiale, appelée aussi approximation de Gauss, consiste à considérer des rayons méridiens proches de l'axe et faiblement inclinés par rapport à l'axe. Soit α₁ l'angle du rayon incident par rapport à l'axe. Les composantes du vecteur unitaire sont :

$$\begin{array}{cc}\hfill & u_1=sin{\alpha }_1\simeq {\alpha }_1\hfill \\ \hfill & v_1=0\hfill \\ \hfill & w_1=cos{\alpha }_1\simeq 1\hfill \end{array}$$

où les expressions approchées sont obtenues en négigeant les termes d'ordre 2.

En négligeant également les termes d'ordre 2 (x₁u₁ et x₁²) dans les expressions de e et f, on obtient :

$$t\simeq z_S-z_1$$

Compte tenu de ces formes approchées, la translation de P₁ à P₂ s'exprime sous la forme matricielle suivante :

$$\left(\begin{array}{c}\hfill n_1{\alpha }_1\hfill \\ \hfill x_2\hfill \end{array}\right)=\left(\begin{array}{cc}\hfill 1\hfill & \hfill 0\hfill \\ \hfill \frac{z_S-z_1}{n_1}\hfill & \hfill 1\hfill \end{array}\right)\left(\begin{array}{c}\hfill n_1{\alpha }_1\hfill \\ \hfill x_1\hfill \end{array}\right)$$

La matrice 2x2 ainsi définie est la matrice de translation entre les plans d'abscisses z₁ et z_S.

La réfraction (ou la réflexion) s'exprime par la relation matricielle suivante :

$$\left(\begin{array}{c}\hfill n_2{\alpha }_2\hfill \\ \hfill x_2\hfill \end{array}\right)=\left(\begin{array}{cc}\hfill 1\hfill & \hfill -C(n_2-n_1)\hfill \\ \hfill 0\hfill & \hfill 1\hfill \end{array}\right)\left(\begin{array}{c}\hfill n_1{\alpha }_1\hfill \\ \hfill x_2\hfill \end{array}\right)$$

On considère à présent un système centré complet formé de N dioptres. Soit z_e la position du dioptre d'entrée et z_s celle du dioptre de sortie. La relation entre ces deux plans s'écrit :

$$\left(\begin{array}{c}\hfill n_s{\alpha }_s\hfill \\ \hfill x_s\hfill \end{array}\right)=\left(\begin{array}{cc}\hfill a_{11}\hfill & \hfill a_{12}\hfill \\ \hfill a_{21}\hfill & \hfill a_{22}\hfill \end{array}\right)\left(\begin{array}{c}\hfill n_e{\alpha }_e\hfill \\ \hfill x_e\hfill \end{array}\right)$$

La matrice de transfert du système A=(a_ij) ne dépend que des dioptres et de leurs positions relatives. On l'obtient en multipliant les matrices de translation et de réfraction définies ci-dessus. Le déterminant de la matrice de transfert est égal à 1.

Pour établir les conditions de stigmatisme, on considère un objet dont la position est définie par :

$$z_o=z_e+t_o$$

et on recherche son image à la position

$$z_i=z_s+t_i$$

L'indice du milieu objet (le plus souvent 1) est noté n_o. Celui du milieu image est noté n_i. Dans le cas où la lumière sort du système dans le sens négatif, ce dernier indice doit être négatif, le plus souvent égal à -1. La même remarque vaut pour la lumière incidente.

Le passage du plan objet au plan image s'écrit :

$$\left(\begin{array}{c}\hfill n_i{\alpha }_i\hfill \\ \hfill x_i\hfill \end{array}\right)=\left(\begin{array}{cc}\hfill b_{11}\hfill & \hfill b_{12}\hfill \\ \hfill b_{21}\hfill & \hfill b_{22}\hfill \end{array}\right)\left(\begin{array}{c}\hfill n_o{\alpha }_o\hfill \\ \hfill x_o\hfill \end{array}\right)$$

où la matrice B est définie par :

$$B=\left(\begin{array}{cc}\hfill 1\hfill & \hfill 0\hfill \\ \hfill \frac{t_i}{n_i}\hfill & \hfill 1\hfill \end{array}\right)\left(\begin{array}{cc}\hfill a_{11}\hfill & \hfill a_{12}\hfill \\ \hfill a_{21}\hfill & \hfill a_{22}\hfill \end{array}\right)\left(\begin{array}{cc}\hfill 1\hfill & \hfill 0\hfill \\ \hfill -\frac{t_o}{n_o}\hfill & \hfill 1\hfill \end{array}\right)$$

Le stigmatisme est obtenu lorsque x_i ne dépend pas de l'angle α_o. Cette condition est réalisée lorsque b₂₁=0, ce qui donne les relations de conjugaison :

$$\frac{t_i}{n_i}=\frac{-a_{21}+a_{22}\frac{t_o}{n_o}}{a_{11}-a_{12}\frac{t_o}{n_o}},\frac{t_o}{n_o}=\frac{-a_{21}-a_{11}\frac{t_i}{n_i}}{-a_{22}-a_{12}\frac{t_i}{n_i}}$$

Le grandissement est :

$$\gamma =\frac{x_i}{x_o}=\frac{1}{a_{11}-a_{12}\frac{t_o}{n_o}}$$

Le foyer objet F_o est le point de l'axe optique dont l'image est à l'infini. Le foyer image F_i est l'image d'un objet situé à l'infini sur l'axe. Leur positions sont données par :

$$t_{Fo}=\frac{n_0{a}_{11}}{a_{12}},t_{Fi}=-\frac{n_i{a}_{22}}{a_{12}}$$

Les plans principaux objet et image sont les plans conjugués donnant un grandissement de 1. Leurs positions sont données par :

$$t_{Ho}=\frac{n_o(a_{11}-1)}{a_{12}},t_{Hi}=\frac{n_i(1-a_{22})}{a_{12}}$$

Enfin les distance focales sont les abscisses des foyers par rapport aux plans principaux correspondants :

$$f_o=t_{Fo}-t_{Ho}=\frac{n_o}{a_{12}},f_i=t_{Fi}-t_{Hi}=-\frac{n_i}{a_{12}}$$

Lorsque l'objet ponctuel est à l'infini, sa position est donnée par une angle α₀. Si le système n'est pas afocal, la position de l'image dans le plan focal est donnée par :

$$x_i=(-\frac{a_{22}}{a_{12}}{a}_{11}+a_{21})n_0{\alpha }_0$$

Dans le cas des télescopes, cette relation permet de calculer la focale effective.

3.a. Translation et réfraction

On considère une surface conique de révolution. L'intersection avec le rayon incident est obtenue à partir de l'équation suivante de la surface conique ([1]) :

$$\frac{C}{2}(x_2^2+y_2^2+S(z_2-z_S{)}^2)-(z_2-z_S)=0$$

C est la courbure au sommet et S le facteur de forme qui s'exprime en fonction de l'excentricité de la conique :

$$S=e^2-1$$

Par exemple, une surface parabolique est obtenue avec S = 0. La surface hyperbolique est obtenue pour S >0 . La surface elliptique est obtenue pour $-1\le S<0$. Pour une hyperbole ou une ellipse, le centre de symétrie de la surface a une abscisse :

$$z_a=-\frac{1}{SC}$$

L'équation à résoudre est :

$$a{t}^2+bt+d=0$$

avec les coefficients suivants :

$$\begin{array}{cc}\hfill & a=\frac{C}{2}(u_1^2+v_1^2+S{w}_1^2)\hfill \\ \hfill & b=C(x_1{u}_1+y_1{v}_1+S(z_1-z_S)w_1)-w_1\hfill \\ \hfill & d=\frac{C}{2}(x_1^2+y_1^2+S{(z_1-z_S)}^2)-(z_1-z_S)\hfill \end{array}$$

La distance du point d'intersection est alors donnée par :

$$t=\frac{-b\pm \sqrt{b^2-4ad}}{2a}$$

Le cas a=0 (qui inclut le cas du dioptre plan) se traite à part.

Pour la parabole, il suffit de retenir la solution positive. Pour une hyperbole ou une ellipse, il faut en plus que l'intersection soit du côté du sommet. Cette dernière condition s'écrit :

$$(z_2-z_a)(z_S-z_a)>0$$

Si deux solutions positives vérifient cette condition, il faut choisir la plus petite.

La cas particulier du dioptre sphérique se retrouve avec S = -1.

La normale unitaire au point d'incidence se calcule en appliquant l'opérateur gradient à l'équation de la conique ([1]). L'angle d'incidence est donné par :

$$cos{\theta }_1=\frac{-C{u}_1{x}_2-C{v}_1{y}_2+w_1(1-CS(z_s-z_2))}{\sqrt{C^2(x_2^2+y_2^2)+{(1-CS(z_s-z_2))}^2}}$$

Si on note D le dénominateur de l'expression précédente, le vecteur unitaire du rayon réfracté est donné par :

$$\begin{array}{cc}\hfill & K=-n_{1}cos{\theta }_1+n_{2}cos{\theta }_2\hfill \\ \hfill & n_2{u}_2=n_1{u}_1-\frac{CK}{D}{x}_2\hfill \\ \hfill & n_2{v}_2=n_1{v}_1-\frac{CK}{D}{y}_2\hfill \\ \hfill & n_2{w}_2=n_1{w}_1+\frac{K}{D}(1-CS(z_s-z_2))\hfill \end{array}$$

4.a. Introduction

Le module Dioptres.py décrit ci-dessous comporte des classes permettant d'effectuer les calculs pour un système centré défini par une liste de dioptres.

4.b. Rayon lumineux

La classe suivante représente un rayon lumineux élémentaire, déterminé par un point et un vecteur unitaire. Les composantes d'un point ou d'un vecteur sont représentées par une liste simple. Le vecteur donné en argument est normalisé.

from math import *
from pylab import *
from IndiceVerre import *

class Rayon:
    def __init__(self,P,U):
        self.P = P
        N = sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2])
        self.U = [U[0]/N,U[1]/N,U[2]/N]
        self.actif = True
        self.longueur = 0.0
    def getPoint(self,t):
        return [self.P[0]+t*self.U[0],self.P[1]+t*self.U[1],self.P[2]+t*self.U[2]]
                

La classe suivante représente un rayon complet, formé d'une suite de rayons élémentaires.

class RayonComplet:
    def __init__(self):
        self.listeR = []
        self.n = 0
    def ajouter(self,R):
        self.listeR.append(R)
        self.n += 1
    def getRayon(self):
        if self.n==0:
            return 0
        else:
            return self.listeR[self.n-1]
            

La méthode suivante fournit le tracé du rayon sous forme de trois listes de coordonnées. L'argument tf précise la longueur de prolongement du dernier rayon élémentaire, s'il est actif.

    def getXYZ(self,tf):
        x=[]
        y=[]
        z=[]
        for k in range(self.n):
           x.append(self.listeR[k].P[0])
           y.append(self.listeR[k].P[1])
           z.append(self.listeR[k].P[2])
        if self.listeR[self.n-1].actif:
            P = self.listeR[self.n-1].getPoint(tf)
            x.append(P[0])
            y.append(P[1])
            z.append(P[2])
        return [x,y,z] 
                

La méthode suivante calcule la longueur optique (longueur multipliée par l'indice) du rayon, en s'arrêtant au rang last avant le dernier rayon.

    def longueurOptique(self,last):
        l = 0.0
        for k in range(self.n-1-last):
            l += self.listeR[k].longueur
        return l
                

4.c. Faisceaux de rayons

Un faisceau est un ensemble de rayons, c'est-à-dire une liste d'instances de la classe RayonComplet. La classe suivante est la classe parente de tous les types de faisceaux.

class Faisceau:
    def __init__(self):
        self.listeRayons = []
        self.n = 0
    def ajouter(self,Rc):
        self.listeRayons.append(Rc)
        self.n += 1
                

La méthode suivante calcule le barycentre des points des derniers rayons actifs du faisceau, qui seront des points d'intersection avec un plan image. Cela permet de localiser le centre de l'image.

    def barycentre(self):
        x = 0
        y = 0
        na = 0
        for Rc in self.listeRayons:
            R = Rc.getRayon()
            if R.actif:
                na +=1 
                x += R.P[0]
                y += R.P[1]
        return [x/na,y/na]
                

Le premier type de faisceau est constitué de N rayons parallèles méridiens, c'est-à-dire contenus dans le plan Oxz, faisant un angle α avec l'axe optique. On doit aussi préciser l'abscisse z de l'intersection du rayon central du faisceau avec l'axe optique et le rayon du faisceau.

class FaisceauParalleleMeridien(Faisceau):
    def __init__(self,nr,alpha,z,r):
        Faisceau.__init__(self)
        dr = 2.0*r/(nr-1)
        cosinus = cos(alpha)
        sinus = sin(alpha)
        for k in range(nr):
            r1 = -r+k*dr
            x = r1*cosinus
            z = -r1*sinus + z
            Rc = RayonComplet()
            Rc.ajouter(Rayon([x,0,z],[sinus,0,cosinus]))
            self.ajouter(Rc)
                

Le type de faisceau suivant est constitué de rayons parallèles répartis sur un cylindre de révolution faisant un angle α par rapport à l'axe optique. On doit préciser l'intersection de l'axe du cylindre avec l'axe optique et le rayon du cylindre. La rotation se fait toujours autours de l'axe y. Remarquer que seuls deux rayons de ce faisceau sont méridiens (si l'angle est non nul), précisément situés dans le plan Oxz. Ces deux rayons sont accessibles par deux méthodes.

class FaisceauParalleleCylindre(Faisceau):
    def __init__(self,nr,alpha,z,r):
        Faisceau.__init__(self)
        N = int(nr/2)
        cosinus = cos(alpha)
        sinus = sin(alpha)
        dtheta = pi/(N-1)
        for k in range(2*N):
            theta = dtheta*k
            xx = r*cos(theta)
            yy = r*sin(theta)
            xp = xx*cosinus
            yp = yy
            zp = -xx*sinus + z
            Rc = RayonComplet()
            Rc.ajouter(Rayon([xp,yp,zp],[sinus,0,cosinus]))
            self.ajouter(Rc)
        

Le faisceau suivant est constitué de rayons parallèles répartis à l'intérieur d'un cylindre de révolution (faisceau collimaté cylindrique). Les rayons se répartissent uniformément sur un plan perpendiculaire au faisceau, en suivant un maillage carré. On précise en argument le nombre de rayons sur un diamètre, l'intersection avec l'axe et le rayon du cylindre.

class FaisceauParallele(Faisceau):
    def __init__(self,nr,alpha,z,r):
        Faisceau.__init__(self)
        dx = 2*r/nr
        r2 = r*r
        cosinus = cos(alpha)
        sinus = sin(alpha)
        for i in range(nr):
            xx = -r+dx*i
            for j in range(nr):
                yy = -r+dx*j
                rr2 = xx*xx+yy*yy
                if (rr2 < r2):
                    xp = xx*cosinus
                    yp = yy
                    zp = -xx*sinus + z
                    Rc = RayonComplet()
                    Rc.ajouter(Rayon([xp,yp,zp],[sinus,0,cosinus]))
                    self.ajouter(Rc)
        

Le faisceau suivant est constitué de rayons méridiens issus d'un point objet dont la position est donnée par (xo,0,zo). Ces rayons éclairent uniformément une pupille de rayon r située en z. Si cette pupille est à gauche de l'objet, celui-ci est virtuel et le faisceau est convergent.

class FaisceauDivergentMeridien(Faisceau):
    def __init__(self,nr,zo,xo,z,r):
        dr = 2*r/(nr-1)
        for k in range(nr):
            r1 = -r+k*dr
            Rc = RayonComplet()
            if z > zo:
                Rc.ajouter(Rayon([xo,0,z0],[r1-xo,0,z-zo]))
            else:
                u = xo-r1
                w = zo-z
                f = 2*w/sqrt(u*u+w*w)
                Rc.ajouter(Rayon([xo-u*f,0,zo-w*f],[u,0,w]))
            self.ajouter(Rc)
            

Le faisceau suivant et analogue au précédent mais éclaire la totalité de la pupille de manière uniforme, par un maillage carré.

class FaisceauDivergent(Faisceau):
    def __init__(self,nr,zo,xo,z,r):
        dx = 2*r/nr;
        r2 = r*r
        for i in range(nr):
            xx = -r+dx*i
            for j in range(nr):
                yy = -r+dx*j
                rr2 = xx*xx+yy*yy
                if rr2 < r2:
                    Rc = RayonComplet()
                    if z > zo:
                        Rc.ajouter(Rayon([xo,0,z0],[xx-xo,yy,z-zo]))
                    else:
                        u = xo-xx
                        v = -yy
                        w = zo-z
                        f = 2*w/sqrt(u*u+v*v+w*w)
                        Rc.ajouter(Rayon([xo-u*f,-v*f,zo-w*f],[u,v,w]))
                self.ajouter(Rc)
            

4.d. Dioptre

La classe suivante est la classe parente des différents types de dioptres (sphériques, coniques, etc). Le dioptre est défini par sa position, sa courbure sur l'axe, ses deux indices et son rayon d'ouverture. Les indices sont des instances de la classe Verre. Le premier indice est celui du milieu traversé par la lumière avant le dioptre (tenir compte du sens de la lumière pour l'attribuer). En revanche la courbure d'un dioptre ne dépend pas du sens de la lumière. Par exemple, un dioptre de courbure positive sera convexe pour une lumière dans le sens direct et concave pour une lumière de sens inverse.

Pour un miroir, l'indice doit changer de signe.

class Dioptre:
    def __init__(self,z,c,n1,n2,r):
        self.z = z
        self.c = c
        self.n1 = n1
        self.n2 = n2
        self.r = r
        self.r2 = r*r
                

Les méthodes suivantes fournissent respectivement la matrice de réfraction pour l'approximation paraxiale et la matrice de translation à partir d'une abscisse z:

    def matriceRefraction(self,raie):
        return [[1,self.c*(self.n1.indice(raie)-self.n2.indice(raie))],[0,1]]
    def matriceTranslation(self,z,raie):
        return [[1,0],[(self.z-z)/self.n1.indice(raie),1]]
            

La méthode suivante calcule le décalage du bord du dioptre par rapport à son sommet, en supposant une forme sphérique :

    def dzBord(self):
        if self.c==0:
            return 0.0
        rc = abs(1/self.c)
        dz = rc-sqrt(rc*rc-self.r*self.r)
        if self.c < 0:
            dz = -dz
        return dz
            

4.e. Dioptre sphérique

La classe suivante représente un dioptre sphérique (ou plan) et effectue les calculs d'intersection avec un rayon et de réfraction. La méthode refraction effectue le calcul de réfraction d'un rayon et renvoit le nouveau rayon. La longueur d'onde en nanomètres est donnée en argument. Si la rayon ne parvient pas au dioptre, ou s'il n'est pas réfracté, il est renvoyé avec la propriété actif = False.

On remarquera que les indices utilisés pour le calcul dépendent du sens de la lumière incidente. Si le dioptre est un miroir, le rayon incident a un indice sont le signe correspond à son sens de propagation. Bien sûr, cela est aussi valable pour le rayon réfléchi.

class Spherique(Dioptre):
    def refraction(self,R1,raie):
        n1 = self.n1.indice(raie)
        n2 = self.n2.indice(raie)
        z1s = R1.P[2]-self.z
        R1.actif = False
        if self.c==0:
            if (R1.U[2]!=0):
                t = -z1s/R1.U[2]
                P2 = R1.getPoint(t)
            else:
                return R1
        else:
            B = (R1.P[0]*R1.U[0]+R1.P[1]*R1.U[1]+z1s*R1.U[2])-R1.U[2]/self.c
            D = (R1.P[0]*R1.P[0]+R1.P[1]*R1.P[1]+z1s*z1s)-2*z1s/self.c
            delta = B*B-D
            if (delta <=0):
                return R1
            sqd = sqrt(delta)
            ta = -B+sqd
            tb = -B-sqd
            zc = self.z+1.0/self.c
            P2 = R1.getPoint(ta)
            a = (P2[2]-zc)*(self.z-zc)
            P2 = R1.getPoint(tb)
            b = (P2[2]-zc)*(self.z-zc)
            if ((a<0)|(ta<0)):
                t = tb
                P2 = R1.getPoint(t)
            elif ((b<0)|(tb<0)):
                t = ta
                P2 = R1.getPoint(t)
            else:
                if (ta < tb):
                    t = ta
                    P2 = R1.getPoint(t)
                else:
                    t = tb
                    P2 = R1.getPoint(t)
        R1.longueur = t*abs(n1)
        if (P2[0]*P2[0]+P2[1]*P2[1] > self.r2):
            return R1
        cosTheta1 = R1.U[2]*(1+self.c*(self.z-P2[2]))-self.c*(R1.U[0]*P2[0]+R1.U[1]*P2[1])
        x = 1-pow(n1/n2,2)*(1-cosTheta1*cosTheta1)
        if x<0:
            return R1
        cosTheta2 = sqrt(x)
        if R1.U[2]>0:
            K = n2*cosTheta2-n1*cosTheta1
        else:
            cosTheta1 = -cosTheta1
            K = -n2*cosTheta2+n1*cosTheta1
        U2 = [0,0,0]
        U2[0] = (n1*R1.U[0]-self.c*K*P2[0])/n2
        U2[1] = (n1*R1.U[1]-self.c*K*P2[1])/n2
        U2[2] = (n1*R1.U[2]-self.c*K*(P2[2]-self.z)+K)/n2
        if n2*n1<0:
            U2[0] = -U2[0]
            U2[1] = -U2[1]
            U2[2] = -U2[2]
        R1.actif = True
        return Rayon(P2,U2)
            

La méthode suivant fournit le profil méridien du dioptre (dans le plan zx) sous forme de deux listes de coordonnées. Le nombre de points est donné en argument.

    def profil(self,npts):
        x = []
        z = []
        if self.c==0:
            x = [-self.r,self.r]
            z = [self.z,self.z]
        else:
            Rc = 1.0/self.c
            zc = self.z + Rc
            alpha = asin(fabs(self.r/Rc))
            if Rc<0:
                alpha1 = -alpha
                alpha2 = alpha
            else:
                alpha1 = pi-alpha
                alpha2 = pi+alpha
            da = (alpha2-alpha1)/(npts-1)
            Rc = fabs(Rc)
            for k in range(npts):
                a = alpha1 + da*k
                z.append(zc+Rc*cos(a))
                x.append(Rc*sin(a))
        return [x,z]
            

4.f. Dioptre conique

Pour ce type de dioptre, il faut préciser le facteur de forme S. On évitera le cas de la courbure nulle qui sera traité avec le dioptre sphérique.

class Conique(Dioptre):
    def __init__(self,z,c,S,n1,n2,r):
        self.S = S
        Dioptre.__init__(self,z,c,n1,n2,r)
    def refraction(self,R1,raie):
        n1 = self.n1.indice(raie)
        n2 = self.n2.indice(raie)
        R1.actif = False
        z1s = R1.P[2]-self.z
        A = self.c/2.0*(R1.U[0]*R1.U[0]+R1.U[1]*R1.U[1]+self.S*R1.U[2]*R1.U[2])
        B = self.c*(R1.P[0]*R1.U[0]+R1.P[1]*R1.U[1]+self.S*z1s*R1.U[2])-R1.U[2]
        D = self.c/2.0*(R1.P[0]*R1.P[0]+R1.P[1]*R1.P[1]+self.S*z1s*z1s)-z1s
        if A==0.0:
            if B==0.0:
                return R1
            t = -D/B
            P2 = R1.getPoint(t)
        else:
            delta = B*B-4*A*D
            if delta<=0.0:
                return R1
            sqd = sqrt(delta)
            ta = (-B+sqd)/(2*A)
            tb = (-B-sqd)/(2*A)
            Pa = R1.getPoint(ta)
            Pb = R1.getPoint(tb)
            if self.S==0.0:
                a=1
                b=1
            else:
                zc = z-1.0/(self.c*self.S)
                a = (Pa[2]-zc)*(self.z-zc)
                b = (Pb[2]-zc)*(self.z-zc)
            if ((a<0)|(ta<0)):
                t = tb
                P2 = R1.getPoint(t)
            elif ((b<0)|(tb<0)):
                t = ta
                P2 = R1.getPoint(t)
            else:
                if (ta < tb):
                    t = ta
                    P2 = R1.getPoint(t)
                else:
                    t = tb
                    P2 = R1.getPoint(t)
        R1.longueur = t*abs(n1)
        r2 = P2[0]*P2[0]+P2[1]*P2[1]
        if ( r2 > self.r2):
            return R1
        az = 1-self.c*self.S*(self.z-P2[2])
        den = sqrt(self.c*self.c*r2+az*az)
        cosTheta1 = (R1.U[2]*az-self.c*(R1.U[0]*P2[0]+R1.U[1]*P2[1]))/den
        x = 1-pow(n1/n2,2)*(1-cosTheta1*cosTheta1)
        if x<0:
            return R1
        cosTheta2 = sqrt(x)
        if R1.U[2]>0:
            K = n2*cosTheta2-n1*cosTheta1
        else:
            cosTheta1 = -cosTheta1
            K = -n2*cosTheta2+n1*cosTheta1
        U2 = [0,0,0]
        U2[0] = (n1*R1.U[0]-self.c*K/den*P2[0])/n2
        U2[1] = (n1*R1.U[1]-self.c*K/den*P2[1])/n2
        U2[2] = (n1*R1.U[2]+K/den*az)/n2
        if n2*n1<0:
            U2[0] = -U2[0]
            U2[1] = -U2[1]
            U2[2] = -U2[2]
        R1.actif = True
        return Rayon(P2,U2)
                

La méthode suivant fournit le profil méridien du dioptre (dans le plan zx) sous forme de deux listes de coordonnées. Le nombre de points est donné en argument.

    def profil(self,npts):
        x = []
        z = []
        if self.c==0:
            x = [-self.r,self.r]
            z = [self.z,self.z]
        else:
            Rc = 1.0/self.c
            zc = self.z + Rc
            alpha = asin(fabs(self.r/Rc))
            if Rc<0:
                alpha1 = -alpha
                alpha2 = alpha
            else:
                alpha1 = pi-alpha
                alpha2 = pi+alpha
            da = (alpha2-alpha1)/(npts-1)
            Rc = fabs(Rc)
            for k in range(npts):
                a = alpha1 + da*k
                z.append(zc+Rc*cos(a))
                x.append(Rc*sin(a))
        return [x,z]
            

4.g. Écran plan

L'écran plan est un disque plan qui intercepte les rayons. Il permet de modéliser une osbtruction centrale. Il est défini par sa position et son rayon.

class Ecran(Dioptre):
    def __init__(self,z,r):
        Dioptre.__init__(self,z,0,IndiceVide(),IndiceVide(),r)
    def refraction(self,R1,raie):
        t = (self.z-R1.P[2])/R1.U[2]
        if (t<0):
            R1.actif = False
            return R1
        na = self.n1.indice(raie)
        R1.longueur = t*abs(na)
        P2 = R1.getPoint(t)
        U2 = R1.U
        R2 = Rayon(P2,U2)
        if (P2[0]*P2[0]+P2[1]*P2[1] < self.r2):
            R2.actif = False
        return R2
    def profil(self,npts):
        x = [-self.r,self.r]
        z = [self.z,self.z]
        return [x,z]
                

La classe suivante effectue la même chose que l'écran met laisse les rayons actifs. Cela permet de l'utiliser pour déterminer l'intersection des rayons avec un plan image.

class PlanImage(Dioptre):
    def __init__(self,z,r):
        Dioptre.__init__(self,z,0,IndiceVide(),IndiceVide(),r)
    def refraction(self,R1,raie):
        t = (self.z-R1.P[2])/R1.U[2]
        na = self.n1.indice(raie)
        R1.longueur = t*abs(na)
        P2 = R1.getPoint(t)
        U2 = R1.U
        R2 = Rayon(P2,U2)
        return R2
    def profil(self,npts):
        x = [-self.r,self.r]
        z = [self.z,self.z]
        return [x,z]
                

4.h. Diaphragme

Le diaphragme intercepte les rayons en dehors d'un disque.

class Diaphragme(Dioptre):
    def __init__(self,z,rmax):
        self.rmax = rmax
        Dioptre.__init__(self,z,0,IndiceVide(),IndiceVide(),r)
    def refraction(self,R1,raie):
        t = (self.z-R1.P[2])/self.U[2]
        if (t<0):
            R1.actif = False
            return R1
        na = self.n1.indice(raie)
        R1.longueur = t*abs(na)
        P2 = R1.getPoint(t)
        U2 = R1.U
        R2 = Rayon(P2,U2)
        if (P2[0]*P2[0]+P2[1]*P2[1] > self.r2):
            R2.actif = False
        return R2
    def profil(self,npts):
        x = [[self.r,self.rmax],[-self.r,-self.rmax]]
        z = [[self.z,self.z],[self.z,self.z]]
        return [x,z]
                

4.i. Intersection avec une sphère

La classe suivante représente un dioptre virtuel (sans effet sur les rayons) qui permet de déterminer l'intersection des rayons avec une sphère, généralement pour déterminer le déphasage avec une surface sphérique de référence. Son comportement est analogue à celui de la classe PlanImage ci-dessus. Le constructeur prend en arguments l'abscisse z du point de la sphère le plus à gauche et les coordonnées PC du centre de la sphère.

class SphereImage(Dioptre):
    def __init__(self,z,PC):
        self.PC = PC
        Dioptre.__init__(self,z,0,IndiceVide(),IndiceVide(),0)
    def refraction(self,R1,raie):
        b = R1.U[0]*(R1.P[0]-self.PC[0])+R1.U[1]*(R1.P[1]-self.PC[1])+R1.U[2]*(R1.P[2]-self.PC[2])
        a0 = R1.P[0]-self.PC[0]
        a1 = R1.P[1]-self.PC[1]
        a2 = R1.P[2]-self.PC[2]
        R = self.PC[2]-self.z
        d = a0*a0+a1*a1+a2*a2-R*R
        delta = b*b-d
        if delta <= 0:
            R1.actif = False
            return R1
        rdelta = sqrt(delta)
        t = min(-b+rdelta,-b-rdelta)
        if t < 0:
            R1.actif = False
            return R1
        na = self.n1.indice(raie)
        R1.longueur = t*abs(na)
        P2 = R1.getPoint(t)
        U2 = R1.U
        R2 = Rayon(P2,U2)
        return R2
    def profil(self,npts):
        return [[0],[0]]
                

4.j. Position sur l'axe optique

La classe suivante représente la position d'un objet ponctuel.

class Position:
    def __init__(self,z):
        self.z = z
        self.infini = False
                

4.k. Système centré

Le système centré complet est constitué d'une liste de dioptres à traiter dans l'ordre. La classe suivante représente ce système et permet d'effectuer les calculs de rayons et d'approximation paraxiale. Par les indices d'entrée et de sortie sont 1.

class SystemeCentre:
    def __init__(self):
        self.listeD = []
        self.nbD = 0
        self.M = [[1,0],[0,1]]
        self.no = Vide()
        self.ni = Vide()
    def ajouter(self,dioptre):
        self.listeD.append(dioptre)
        self.nbD += 1
    def ajouterListe(self,liste):
        for d in liste:
            self.ajouter(d)
    def refraction(self,rayonComplet,raie):
        k = 0
        R1 = rayonComplet.getRayon()
        actif = R1.actif
        while k < self.nbD and actif:
            R2 = self.listeD[k].refraction(R1,raie)
            rayonComplet.ajouter(R2)
            R1 = R2
            k += 1
            actif = R2.actif
    def refractionFaisceau(self,faisceau,raie):
        for Rc in faisceau.listeRayons:
            self.refraction(Rc,raie)
                

Les méthodes suivantes effectuent les calculs de l'approximation paraxiale. La première méthode calcule la matrice de transfert du système; elle doit être appelée avant les autres méthodes.

    def matriceTransfert(self,raie):
        self.M = self.listeD[0].matriceRefraction(raie)
        zd = self.listeD[0].z
        k = 1
        while k < self.nbD:
            MT = self.listeD[k].matriceTranslation(zd,raie)
            self.M = self.multiplier(MT,self.M)
            MR = self.listeD[k].matriceRefraction(raie)
            self.M = self.multiplier(MR,self.M)
            zd = self.listeD[k].z
            k += 1
        return self.M
    def multiplier(self,A,B):
        C = [[0,0],[0,0]]
        C[0][0] = A[0][0]*B[0][0]+A[0][1]*B[1][0]
        C[0][1] = A[0][0]*B[0][1]+A[0][1]*B[1][1]
        C[1][0] = A[1][0]*B[0][0]+A[1][1]*B[1][0]
        C[1][1] = A[1][0]*B[0][1]+A[1][1]*B[1][1]
        return C
                

La méthode suivante calcule la position d'une image d'un objet.

    def image(self,PObjet,raie):
        PImage = Position(0)
        if PObjet.infini:
            if self.M[0][1]==0:
                PImage.infini = True
            else:
                ti = -self.ni.indice(raie)*self.M[1][1]/self.M[0][1]
                PImage.z = self.listeD[self.nbD-1].z+ti
        else:
            to = PObjet.z-self.listeD[0].z
            den = self.M[0][0]-self.M[0][1]*to/self.no.indice(raie)
            if den==0:
                PImage.infini = True
            else:
                ti = self.ni.indice(raie)*(-self.M[1][0]+self.M[1][1]*to/self.no.indice(raie))/den
                PImage.z = self.listeD[self.nbD-1].z+ti
        return PImage
                
                

La méthode suivante calcule le grandissement lorsque l'objet et l'image sont à distance finie. Si l'objet est à l'infini et l'image à distance finie, la valeur renvoyée est le rapport de la position de l'image sur l'angle définissant la position de l'objet. Si l'objet et l'image sont à l'infini (système afocal), la valeur renvoyée est le rapport des deux angles, c'est-à-dire le grossissement. Si l'objet est à distance finie et l'image est à l'infini, la valeur renvoyée est le rapport de l'angle donnant la position de l'image sur la position de l'objet.

    def grandissement(self,PObjet,raie):
        if PObjet.infini:
            if self.M[0][1]==0:
                return self.M[0][0]*self.no.indice(raie)/self.ni.indice(raie)
            return (-self.M[1][1]/self.M[0][1]*self.M[0][0]+self.M[1][0])*self.no.indice(raie)
        else:
            to = PObjet.z-self.listeD[0].z
            den = self.M[0][0]-self.M[0][1]*to/self.no.indice(raie)
            if den==0:
                return self.M[0][1]/self.ni.indice(raie)
            else:
                return 1.0/den
                

La méthode suivante calcule la position des foyers objet et image, des plans principaux objet et image.

    def foyers(self,raie):
        Fo = Position(0)
        Fi = Position(0)
        Ho = Position(0)
        Hi = Position(0)
        if self.M[0][1]==0:
            Fo.infini = True
            Fi.infini = True
            Ho.infini = True
            Hi.infini = True
        else:
            Fo.z = self.listeD[0].z + self.no.indice(raie)*self.M[0][0]/self.M[0][1]
            Fi.z = self.listeD[self.nbD-1].z - self.ni.indice(raie)*self.M[1][1]/self.M[0][1]
            Ho.z = self.listeD[0].z + self.no.indice(raie)*(self.M[0][0]-1)/self.M[0][1]
            Hi.z = self.listeD[self.nbD-1].z + self.ni.indice(raie)*(1-self.M[1][1])/self.M[0][1]
        return [Fo,Fi,Ho,Hi]
                

La méthode suivante effectue le tracé des profils des dioptres :

    def tracerProfils(self,npts,color):
        for d in self.listeD:
            [x,z] = d.profil(npts)
            plot(z,x,color='k')
                

La méthode suivante calcule le point de convergence des rayons parallèles à l'axe dont la distance à celui-ci est donnée en argument. Cette fonction permet le calcul des aberrations sphériques et chromatiques longitudinales. Si la distance est nulle, le calcul se fait avec la matrice de l'approximation paraxiale. L'abscisse z₀ fournie en argument correspond au point de départ des rayons, qui doit être à gauche du premier dioptre.

    def pointConvergence(self,z0,r,raie):
        if r==0.0:
            self.matriceTransfert(raie)
            [Fo,Fi,Ho,Hi] = self.foyers(raie)
            return Fi.z
        else:
            R = Rayon([r,0,z0],[0,0,1.0])
            Rc = RayonComplet()
            Rc.ajouter(R)
            self.refraction(Rc,raie)
            R2 = Rc.getRayon()
            t = -R2.P[0]/R2.U[0]
            z = R2.P[2]+R2.U[2]*t
            return z
                

4.l. Exemple : lentille biconvexe

from pylab import *
from Dioptres import *
cat = CatalogueVerre()
n1 = Vide()
n2 = cat.verre['N-BK7']
dioptre1 = Spherique(-0.1,1.0,n1,n2,0.4)
dioptre2 = Spherique(0.1,-1.0,n2,n1,0.4)
sys = SystemeCentre()
sys.ajouter(dioptre1)
sys.ajouter(dioptre2)
figure()
zObjet = -2.0
uxList = [0,0.02,0.04,0.06,0.08,0.1]
for ux in uxList:
    R = Rayon([0,0,zObjet],[ux,0,1.0])
    Rc = RayonComplet()
    Rc.ajouter(R)
    sys.refraction(Rc,'d')
    XYZ = Rc.getXYZ(3)
    plot(XYZ[2],XYZ[0],color='r')
[x,z] = dioptre1.profil(20)
plot(z,x,color='b')
[x,z] = dioptre2.profil(20)
plot(z,x,color='b')
xlabel('z')
ylabel('x')
axis([-3,3,-0.5,0.5])
grid(True)
                

plotA.pdf

PObjet = Position(-2.0)
sys.matriceTransfert('d')
PImage = sys.image(PObjet,'d')
[Fo,Fi,Ho,Hi] = sys.foyers('d')
                

print(PImage.z)
--> 2.071235202871021

print(sys.grandissement(PObjet,'d'))
--> -1.0361920820654342

print([Fo.z,Fi.z,Ho.z,Hi.z])
--> [-1.0333653626064212, 1.0333653626064212, -0.031746205089002985, 0.031746205089002985]

4.m. Exemple : miroir parabolique

n1 = Vide()
n2 = Vide()
n2.negatif()
miroir = Conique(0,-1.0,0.0,n1,n2,0.5)
sys = SystemeCentre()
sys.ajouter(miroir)
figure()
f = FaisceauParalleleMeridien(10,0,-1.0,0.4)
sys.refractionFaisceau(f,'d')
for Rc in f.listeRayons:
    XYZ = Rc.getXYZ(0.6)
    plot(XYZ[2],XYZ[0],color='r')
xlabel('z') 
ylabel('x')
axis([-1,0.5,-0.5,0.5])
grid(True)
                

plotB.pdf

sys.ni.negatif()
sys.matriceTransfert('d')
[Fo,Fi,Ho,Hi] = sys.foyers('d')
                

print([Fo.z,Fi.z,Ho.z,Hi.z])
--> [-0.5, -0.5, 0.0, 0.0]