Table des matières Python

Dioptres centrés : optique géométrique

1. Introduction

Un système optique centré est constitué de lentilles et miroirs dont les surfaces présentent une symétrie de révolution autour d'un axe. On considère ici le calcul des rayons lumineux pour ce type de système dans le cadre de l'optique géométrique. L'approximation paraxiale permet d'obtenir la relation de conjugaison du système. Le calcul complet des rayons permet de mettre en évidence les différentes formes d'aberrations géométriques.

2. Dioptres plans et sphériques

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 n1 et n2. On note S le sommet, c'est-à-dire l'intersection du dioptre avec l'axe optique Oz, et C son centre de courbure.

figureA.svgFigure pleine page

Par convention, l'indice n1 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=1R=1SC¯

Soit un rayon lumineux partant d'un point P1(x1,y1,z1) dans une direction donnée par le vecteur unitaire U1(u1,v1,w1). Ce rayon rencontre la surface du dioptre au point P2(x2,y2,z2) où la normale unitaire est N, avec un angle d'incidence θ1. Le rayon réfracté est défini par le vecteur unitaire U2(u2,v2,w2).

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

x2=x1+tu1y2=y1+tv1z2=z1+tw1

t est la distance parcourue de P1 à P2. 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 z2=zs et la distance est donc:

t=zs-z1w1

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 :

x22+y22+(z2-zS-R)2=R2

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

t2+2bt+d=0

Les coefficients de cette équation sont :

b=x1u1+y1v1+(z1-zS)w1-Rw1d=x12+y12+(z1-zS)2-2R(z1-zS)

Le discriminant de l'équation est :

Δ=b2-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±b2-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 :

(z2-zc)(zS-zc)>0

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

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

N=CP2R

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 :

KN=n1U1-n2U2

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

K=-n1cosθ1+n2cosθ2

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

-U1N=cosθ1=w1R+zS-z2R-u1x2R-v1y2R

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

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

sinθ2=n1n2sinθ1cosθ2=1-sin2θ2

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 :

n2u2=n1u1-CKx2n2v2=n1v1-CKy2n2w2=n1w1-CK(z2-R-zS)=n1w1-CK(z2-zs)+K

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

n2=-n1

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

n1=1n2=-1

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 α1 l'angle du rayon incident par rapport à l'axe. Les composantes du vecteur unitaire sont :

u1=sinα1α1v1=0w1=cosα11

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

En négligeant également les termes d'ordre 2 (x1u1 et x12) dans les expressions de e et f, on obtient :

tzS-z1

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

(n1α1x2)=(10zS-z1n11)(n1α1x1)

La matrice 2x2 ainsi définie est la matrice de translation entre les plans d'abscisses z1 et zS.

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

(n2α2x2)=(1-C(n2-n1)01)(n1α1x2)

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

(nsαsxs)=(a11a12a21a22)(neαexe)

La matrice de transfert du système A=(aij) 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 :

zo=ze+to

et on recherche son image à la position

zi=zs+ti

L'indice du milieu objet (le plus souvent 1) est noté no. Celui du milieu image est noté ni. 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 :

(niαixi)=(b11b12b21b22)(noαoxo)

où la matrice B est définie par :

B=(10tini1)(a11a12a21a22)(10-tono1)

Le stigmatisme est obtenu lorsque xi ne dépend pas de l'angle αo. Cette condition est réalisée lorsque b21=0, ce qui donne les relations de conjugaison :

tini=-a21+a22tonoa11-a12tono,tono=-a21-a11tini-a22-a12tini

Le grandissement est :

γ=xixo=1a11-a12tono

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

tFo=n0a11a12,tFi=-nia22a12

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

tHo=no(a11-1)a12,tHi=ni(1-a22)a12

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

fo=tFo-tHo=noa12,fi=tFi-tHi=-nia12

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

xi=(-a22a12a11+a21)n0α0

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

3. Dioptres asphériques coniques

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]) :

C2(x22+y22+S(z2-zS)2)-(z2-zS)=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=e2-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 -1S<0 . Pour une hyperbole ou une ellipse, le centre de symétrie de la surface a une abscisse :

za=-1SC

L'équation à résoudre est :

at2+bt+d=0

avec les coefficients suivants :

a=C2(u12+v12+Sw12)b=C(x1u1+y1v1+S(z1-zS)w1)-w1d=C2(x12+y12+S(z1-zS)2)-(z1-zS)

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

t=-b±b2-4ad2a

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 :

(z2-za)(zS-za)>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θ1=-Cu1x2-Cv1y2+w1(1-CS(zs-z2))C2(x22+y22)+(1-CS(zs-z2))2

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

K=-n1cosθ1+n2cosθ2n2u2=n1u1-CKDx2n2v2=n1v1-CKDy2n2w2=n1w1+KD(1-CS(zs-z2))

4. Module python

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 z0 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)
                
plotAplotA.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)
                
plotBplotB.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]
Références
[1]  A. Nussbaum,  Optical system design,  (Prentice Hall PTR, 1998)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.