Accueil > Équipes > S2I > Eric Anterrieu > Ouvrages > Simulations numériques en optique

 

Simulations numériques en optique

Par Eric Anterrieu - 1er/12/2005

 

AVIS IMPORTANT
Les programmes présentés ici sont destinés à illustrer les simulations numériques de l’annexe 3 du livre :
Optique : fondements et applications
de José-Philippe PÉREZ avec la collaboration de Éric ANTERRIEU, DUNOD éditeur (réf. ISBN : 2-10-048497-4).

L’utilisation de ces programmes, écrits pour le logiciel de calcul scientifique MATLAB, est libre de droit et n’engage nullement la responsabilité des auteurs ni celle de l’éditeur ; elle ne peut en aucun cas se faire à des fins commerciales sans l’accord des auteurs et de l’éditeur, ce qui n’exclut donc pas d’utiliser ces programmes en classe ou en séance de travaux pratiques ou dirigés, à la condition expresse que la licence MATLAB soit à jour.

I.- Structures de données et fonctions de base

On utilise les capacités de programmation orienté objet (POO) de MATLAB afin de mettre en oeuvre des classes qui renferment toutes les caractéristiques physiques (attributs ou propriétés en POO) et toutes les actions que l’on peut effectuer (fonctions membre ou messages) sur des systèmes simples comme le dioptre, le miroir et le diaphragme. On associe des objets de ces trois classes en utilisant la fonction sysopt.m comme on associe des dioptres, des miroirs et des diaphragmes sur un banc optique :


function obj=sysopt(varargin)
% obj=SYSOPT(varargin) returns a cell array where each cell contains the optical
% characteristics of variable number of optical subsystem.
% Subsystems should be passed in the order they are hitted by the light.
%  » varargin = optical subsystems (variable number)
%  « obj = optical system
%
% © Eric Anterrieu (2004)
%
% See also: DIOPTRE/DIOPTRE, MIROR/MIROR, DIAPHRAGM/DIAPHRAGM.
%
if (nargin < 1), error('Wrong number of input arguments.'); end
%
for k=1:length(varargin)
  obj{k+1}=varargin{k};
end
obj{1}=[get(varargin{1},'FrontIndex') ...
        get(varargin{end},'BackIndex') ...
        get(varargin{1},'Position') ...
        get(varargin{end},'Position')];
%
return;
sysopt.m

I.1.- Éléments constitutifs

a) La classe dioptre

Le fichier dioptre.m dans le dossier @dioptre contient la définition de la classe dioptre :


function obj=dioptre(varargin)
% DIOPTRE is the constructor function for dioptre object.
% The dioptre class is primarily meant to be used to store optical characteristics of a dioptre.
% obj=DIOPTRE(varargin)
% Called with no arguments, the constructor returns a default dioptre object.
% Called with one argument that is a dioptre object, the object is simply returned.
% Called with five arguments (name,Z,R,No,Ni), the constructor returns a new dioptre object.
%  » name = name of the dioptre in the optical system
%  » Z = position of vertex on optical axis [m]
%  » R = radius of curvature [m]
%  » No = refractive index of front medium
%  » Ni = refractive index of back medium
%  « obj = dioptre object
%
% © Eric Anterrieu (2004)
%
% See also: DIOPTRE/DISPLAY, DIOPTRE/GET.
%
switch (nargin)
case 0
  obj.name = '';
  obj.Z = 0;
  obj.R = 0;
  obj.N = [0, 0];
  obj = class(obj,'dioptre');
case 1
  if (isa(varargin{1},'dioptre'))
    obj = varargin{1};
  else
    error('Wrong argument type.');
   end
case 5
  if (varargin{4} < 1), error('Refractive index should be greater than one.');  end
  if (varargin{5} < 1), error('Refractive index should be greater than one.');  end
  obj.name = varargin{1};
  obj.Z    = varargin{2};
  obj.R    = varargin{3};
  obj.N    = [varargin{4}, varargin{5}];
  obj = class(obj,'dioptre');
otherwise
  error('Wrong number of input arguments.');
 end
%
return;
@dioptre/dioptre.m

Dans le dossier @dioptre, on trouve la sur-définition de fonctions pour afficher/modifier les valeurs des propriétés d’un objet de la classe dioptre...


function display(obj)
% DISPLAY(obj) displays the optical characteristics of a dioptre.
%  » obj = dioptre object
%
% © Eric Anterrieu (2004)
%
% See also: DIOPTRE/DIOPTRE, DIOPTRE/GET, DIOPTRE/SET.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
disp(sprintf('               name of the dioptre: %s',  obj.name));
disp(sprintf('position of vertex on optical axis: %g m',obj.Z));
disp(sprintf('               radius of curvature: %g m',obj.R));
disp(sprintf('  refractive index of front medium: %g',  obj.N(1)));
disp(sprintf('   refractive index of back medium: %g',  obj.N(2)));
%
return;
@dioptre/display.m


function obj=set(obj,varargin)
% obj=SET(obj,property,value,...) sets dioptre properties to specified values.
%  » obj = dioptre object
%  » property = dioptre property
%  » val = value of the property
% Valid properties are: 'Name', 'Position', 'Curvature', 'FrontIndex', 'BackIndex'
%
% © Eric Anterrieu (2004)
%
% See also: DIOPTRE/DIOPTRE, DIOPTRE/GET, DIOPTRE/DISPLAY.
%
property_argin = varargin;
while length(property_argin) >= 2
  property = property_argin{1};
  val = property_argin{2};
  property_argin = property_argin(3:end);
  switch property
  case 'Name'
    obj.name = val;
  case 'Position'
    obj.Z = val;
  case 'Curvature'
    obj.R = val;
  case 'FrontIndex'
    obj.N(1) = val;
  case 'BackIndex'
    obj.N(2) = val;
  otherwise
    error([property,' is not a valid dioptre property.']);
  end
end
%
return;
@dioptre/set.m


function val=get(obj,property)
% val=GET(obj,property) returns value of a dioptre property.
%  » obj = dioptre object
%  » property = dioptre property
%  « val = value of the property
% Valid properties are: 'Name', 'Position', 'Curvature', 'FrontIndex', 'BackIndex'
%
% © Eric Anterrieu (2004)
%
% See also: DIOPTRE/DIOPTRE, DIOPTRE/SET, DIOPTRE/DISPLAY.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
ifisa(obj,'dioptre')), error('Wrong argument type.');  end
ifisa(property,'char')), error('Wrong argument type.');  end
%
switch (property)
case 'Name'
  val = obj.name;
case 'Position'
  val = obj.Z;
case 'Curvature'
  val = obj.R;
case 'FrontIndex'
  val = obj.N(1);
case 'BackIndex'
  val = obj.N(2);
otherwise
  error([property,' is not a valid dioptre property.']);
end
%
return;
@dioptre/get.m

Dans le dossier @dioptre, on trouve aussi des fonctions internes à la classe dioptre...


function mat=matrix(obj)
% mat=MATRIX(obj) returns the refraction matrix of a dioptre.
%  » obj = dioptre object
%  « mat = refraction matrix
%
% © Eric Anterrieu (2004)
%
% See also: DIOPTRE/DIOPTRE, DIOPTRE/POWER, DIOPTRE/GET.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
mat      =  eye(2);
mat(2,1) = -power(obj);
%
return;
@dioptre/matrix.m


function V=power(obj)

% V=POWER(obj) returns the refractive power of a dioptre.
%  » obj = dioptre object
%  « V = refractive power of the dioptre
%
% © Eric Anterrieu (2004)
%
% See also: DIOPTRE/DIOPTRE, DIOPTRE/MATRIX, DIOPTRE/GET.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
if (isinf(obj.R))
  V = 0;
else
  V = (obj.N(2)-obj.N(1))/obj.R;
end
%
return;
@dioptre/power.m


function mat=propagation(obj1,obj2)
% mat=PROPAGATION(obj1,obj2) returns the propagation matrix from a dioptre to
% another optical subsystem.
%  » obj1 = dioptre object
%  » obj2 = optical subsystem
%  « mat = propagation matrix
%
% © Eric Anterrieu (2004)
%
% See also: DIOPTRE/DIOPTRE, DIOPTRE/MATRIX, DIOPTRE/GET.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
mat = eye(2);
if (isa(obj2,'dioptre') == 1)
  mat(1,2)=abs(obj2.Z - obj1.Z)/obj1.N(2);
elseif (isa(obj2,'mirror') == 1)
  mat(1,2)=abs(get(obj2,'Position') - obj1.Z)/obj1.N(2);
elseif (isa(obj2,'diaphragm') == 1)
  mat(1,2)=abs(get(obj2,'Position') - obj1.Z)/obj1.N(2);
end
%
return;
@dioptre/propagation.m

Dans le dossier @dioptre, on trouve aussi des fonctions destinées au tracé de rayons...


function h=plot(obj,clipping,varargin)
% h=PLOT(obj,clipping,varargin) draws a dioptre and returns a handle to the line object.
%  » obj = dioptre object
%  » clipping = optical subsystem
%  » varargin = line attributs
%  « h = handle to line object
%
% © Eric Anterrieu (2004)
%
% See also: PLOTSYSOPT, DIOPTRE/DIOPTRE.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
if (isinf(obj.R))
  x=[ 1 1]*obj.Z;
  y=[-1 1]*clipping;
else
  [x,y]=circle(obj.Z+obj.R,obj.R,clipping);
end
if (nargout == 0),   plot(x,y,varargin{:}); end
if (nargout == 1), h=plot(x,y,varargin{:}); end
%
return;

function [x,y]=circle(C,R,clipping) % if (clipping > abs(R)) alpha=90; else alpha=asin(clipping/abs(R))*180/pi; else if (R > 0) a=linspace(180-alpha,180+alpha)*pi/180; else a=linspace(-alpha,alpha)*pi/180; end x=abs(R)*cos(a)+C; y=abs(R)*sin(a); return;
@dioptre/plot.m


function Rray=ray(obj,Iray)
% Rray=RAY(obj,Iray) returns the point hitted by an incidence ray on a dioptre
% as well as the direction of the refracted ray.
%  » obj = dioptre object
%  » Iray(3) = incident ray
%  « Rray(3) = refracted ray
% First and second components of the ray are the co-ordinates of a point on the ray,
% while the third one is the angle made by the ray with the optical axis.
%
% © Eric Anterrieu (2004)
%
% See also: DIOPTRE/DIOPTRE.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
if (isnan(Iray(3))), Rray=[NaN;NaN;NaN]; return; end
%
Z = obj.Z;
R = obj.R;
if (abs(Iray(3)) < pi/2), N1 = obj.N(1); N2 = obj.N(2); % propagation from left to right
else                      N1 = obj.N(2); N2 = obj.N(1); % propagation from right to right
end
% U1
U1 = [cos(Iray(3)), sin(Iray(3))];
% I
a = tan(Iray(3));
b = Iray(2) - a*Iray(1);
if ˜isinf(R)
  [Rray(1),Rray(2)]=Icircle(a,b,Z+R,R);
else
  Rray(1) = Z;
  Rray(2) = a*Rray(1) + b;
end
if isnan(Rray(1)), Rray=[NaN;NaN;NaN]; return; end
% N
if ˜isinf(R)
  [N(1),N(2)]=Ncircle(Rray(1),Rray(2),Z+R);
else
  N = [1, 0];
end
% i1
i1 = acos(dot(U1,N));
if (i1 >= pi/2), i1 = pi-i1; end
% i2
if (sin(i1) >= N2/N1), Rray=[NaN;NaN;NaN]; return; end
i2 = asin(N1*sin(i1)/N2);
if (i2 < 0), i2 = -i2; end
% U2
k = N2*cos(i2) - N1*cos(i1);
U2 = (N1*U1+sign(dot(N,U1))*k*N)/N2;
U2 = U2/norm(U2);
Rray(3) = atan2(U2(2),U2(1));
return;

function [xI,yI]=Icircle(a,b,xC,R) % Return the intersection of a line with a circle. % delta = (1 + a^2)*R^2 - (a*xC + b)^2; if (delta < 0) xI=NaN; yI=NaN; else xI = (xC - a*b - sign(R)*sqrt(delta))/(1 + a^2); yI = a*xI + b; end return;
function [xN,yN]=Ncircle(xI,yI,xC) % Return the normal to a circle. % OC = [xC; 0]; OI = [xI; yI]; IC = OC-OI; IC = IC/norm(IC); xN = IC(1); yN = IC(2); return;
@dioptre/ray.m

b) La classe miroir

Le fichier mirror.m dans le dossier @mirror contient la définition de la classe mirror :


function obj=mirror(varargin)
% MIRROR is the constructor function for mirror object.
% The mirror class is primarily meant to be used to store optical characteristics of a mirror.
% obj=MIRROR(varargin)
% Called with no arguments, the constructor returns a default mirror object.
% Called with one argument that is a mirror object, the object is simply returned.
% Called with four arguments (name,Z,R,N), the constructor returns a new mirror object.
%  » name = name of the dioptre in the optical system
%  » Z = position of vertex on optical axis [m]
%  » R = radius of curvature [m]
%  » N = refractive index of medium
%  « obj = mirror object
%
% © Eric Anterrieu (2004)
%
% See also: MIRROR/DISPLAY, MIRROR/GET.
%
switch (nargin)
case 0
  obj.name = '';
  obj.Z = 0;
  obj.R = 0;
  obj.N = [0, 0];
  obj = class(obj,'mirror');
case 1
  if (isa(varargin{1},'mirror'))
    obj = varargin{1};
  else
    error('Wrong argument type.');
   end
case 4
  if (varargin{4} < 1), error('Refractive index should be greater than one.');  end
  obj.name = varargin{1};
  obj.Z    = varargin{2};
  obj.R    = varargin{3};
  obj.N    = [varargin{4}, varargin{4}];
  obj = class(obj,'mirror');
otherwise
  error('Wrong number of input arguments.');
 end
%
return;
@mirror/mirror.m

Dans le dossier @mirror, on trouve la sur-définition de fonctions pour afficher/modifier les valeurs des propriétés d’un objet de la classe mirror...


function display(obj)
% DISPLAY(obj) displays the optical characteristics of a mirror.
%  » obj = mirror object
%
% © Eric Anterrieu (2004)
%
% See also: MIRROR/MIRROR, MIRROR/GET, MIRROR/SET.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
disp(sprintf('                name of the mirror: %s',  obj.name));
disp(sprintf('position of vertex on optical axis: %g m',obj.Z));
disp(sprintf('               radius of curvature: %g m',obj.R));
disp(sprintf('  refractive index of front medium: %g',  obj.N(1)));
%
return;
@mirror/display.m


function obj=set(obj,varargin)

% obj=SET(obj,property,value,...) sets mirror properties to specified values.
%  » obj = mirror object
%  » property = mirror property
%  » val = value of the property
% Valid properties are: 'Name', 'Position', 'Curvature', 'FrontIndex'
%
% © Eric Anterrieu (2004)
%
% See also: MIRROR/MIRROR, MIRROR/GET, MIRROR/DISPLAY.
%
property_argin = varargin;
while length(property_argin) >= 2
  property = property_argin{1};
  val = property_argin{2};
  property_argin = property_argin(3:end);
  switch property
  case 'Name'
    obj.name = val;
  case 'Position'
    obj.Z = val;
  case 'Curvature'
    obj.R = val;
  case 'FrontIndex'
    obj.N(1) = val;
    obj.N(2) = val;
  otherwise
    error([property,' is not a valid mirror property.']);
  end
end
%
return;
@mirror/set.m


function val=get(obj,property)
% val=GET(obj,property) returns value of a mirror property.
%  » obj = mirror object
%  » property = mirror property
%  « val = value of the property
% Valid properties are: 'Name', 'Position', 'Curvature', 'FrontIndex', 'BackIndex'
%
% © Eric Anterrieu (2004)
%
% See also: MIRROR/MIRROR, MIRROR/SET, MIRROR/DISPLAY.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
ifisa(obj,'mirror')), error('Wrong argument type.');  end
ifisa(property,'char')), error('Wrong argument type.');  end
%
switch (property)
case 'Name'
  val = obj.name;
case 'Position'
  val = obj.Z;
case 'Curvature'
  val = obj.R;
case 'FrontIndex'
  val = obj.N(1);
case 'BackIndex'
  val = obj.N(2);
otherwise
  error([property,' is not a valid mirror property.']);
end
%
return;
@mirror/get.m

Dans le dossier @mirror, on trouve aussi des fonctions internes à la classe mirror...


function mat=matrix(obj,dir)
% mat=MATRIX(obj) returns the reflection matrix of a mirror.
%  » obj = mirror object
%  » dir = direction of propagation of light: '+' or '-'
%  « mat = reflection matrix
%
% © Eric Anterrieu (2004)
%
% See also: MIRROR/MIRROR, MIRROR/POWER, MIRROR/GET.
%
if (nargin < 1), error('Wrong number of input arguments.'); end
if (nargin > 2), error('Wrong number of input arguments.'); end
%
mat = eye(2);
if (nargin == 1)
  mat(2,1) = -power(obj,'+');
else
  mat(2,1) = -power(obj,dir);
end
%
return;
@mirror/matrix.m


function V=power(obj,dir)
% V=POWER(obj) returns the power of a mirror.
%  » obj = mirror object
%  » dir = direction of propagation of light: '+' or '-'
%  « V = power of the mirror
%
% © Eric Anterrieu (2004)
%
% See also: MIRROR/MIRROR, MIRROR/MATRIX, MIRROR/GET.
%
if (nargin < 1), error('Wrong number of input arguments.'); end
if (nargin > 2), error('Wrong number of input arguments.'); end
%
if (isinf(obj.R))
  V = 0;
else
  V = -2*obj.N(1)/obj.R;
  if (nargin == 2)
    if (dir == '-'), V = -V; end
  end
end
%
return;
@mirror/power.m


function mat=propagation(obj1,obj2)
% mat=PROPAGATION(obj1,obj2) returns the propagation matrix from a mirror to
% another optical subsystem.
%  » obj1 = mirror object
%  » obj2 = optical subsystem
%  « mat = propagation matrix
%
% © Eric Anterrieu (2004)
%
% See also: MIRROR/MIRROR, MIRROR/MATRIX, MIRROR/GET.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
mat = eye(2);
if (isa(obj2,'mirror') == 1)
  mat(1,2)=abs(obj2.Z - obj1.Z)/obj1.N(2);
elseif (isa(obj2,'dioptre') == 1)
  mat(1,2)=abs(get(obj2,'Position') - obj1.Z)/obj1.N(2);
elseif (isa(obj2,'diaphragm') == 1)
  mat(1,2)=abs(get(obj2,'Position') - obj1.Z)/obj1.N(2);
end
%
return;
@mirror/propagation.m

Dans le dossier @mirror, on trouve aussi des fonctions destinées au tracé de rayons...


function h=plot(obj,clipping,varargin)
% h=PLOT(obj,clipping,varargin) draws a mirror and returns a handle to the line object.
%  » obj = mirror object
%  » clipping = optical subsystem
%  » varargin = line attributs
%  « h = handle to line object
%
% © Eric Anterrieu (2004)
%
% See also: PLOTSYSOPT, MIRROR/MIRROR.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
if (isinf(obj.R))
  x=[ 1 1]*obj.Z;
  y=[-1 1]*clipping;
else
  [x,y]=circle(obj.Z+obj.R,obj.R,clipping);
end
if (nargout == 0),   plot(x,y,varargin{:}); end
if (nargout == 1), h=plot(x,y,varargin{:}); end
%
return;

function [x,y]=circle(C,R,clipping) % if (clipping > abs(R)) alpha=90; else alpha=asin(clipping/abs(R))*180/pi; else if (R > 0) a=linspace(180-alpha,180+alpha)*pi/180; else a=linspace(-alpha,alpha)*pi/180; end x=abs(R)*cos(a)+C; y=abs(R)*sin(a); return;
@mirror/plot.m


function Rray=ray(obj,Iray)
% Rray=RAY(obj,Iray) returns the point hitted by an incidence ray on a mirror
% as well as the direction of the reflected ray.
%  » obj = mirror object
%  » Iray(3) = incident ray
%  « Rray(3) = reflected ray
% First and second components of the ray are the co-ordinates of a point on the ray,
% while the third one is the angle made by the ray with the optical axis.
%
% © Eric Anterrieu (2004)
%
% See also: MIRROR/MIRROR.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
if (isnan(Iray(3))), Rray=[NaN;NaN;NaN]; return; end
%
Z = obj.Z;
R = obj.R;
% U1
U1 = [cos(Iray(3)), sin(Iray(3))];
% I
a = tan(Iray(3));
b = Iray(2) - a*Iray(1);
if ˜isinf(R)
  [Rray(1),Rray(2)]=Icircle(a,b,Z+R,R);
else
  Rray(1) = Z;
  Rray(2) = a*Rray(1) + b;
end
if isnan(Rray(1)), Rray=[NaN;NaN;NaN]; return; end
% N
if ˜isinf(R)
  [N(1),N(2)]=Ncircle(Rray(1),Rray(2),Z+R);
else
  N = [1, 0];
end
% i1
i1 = acos(dot(U1,N));
if (i1 >= pi/2), i1 = pi-i1; end
% i2
i2 = i1;
% U2
k = -cos(i2) - cos(i1);
U2 = (U1+sign(dot(N,U1))*k*N);
U2 = U2/norm(U2);
Rray(3) = atan2(U2(2),U2(1));
return;

function [xI,yI]=Icircle(a,b,xC,R) % Return the intersection of a line with a circle. % delta = (1 + a^2)*R^2 - (a*xC + b)^2; if (delta < 0) xI=NaN; yI=NaN; else xI = (xC - a*b - sign(R)*sqrt(delta))/(1 + a^2); yI = a*xI + b; end return;
function [xN,yN]=Ncircle(xI,yI,xC) % Return the normal to a circle. % OC = [xC; 0]; OI = [xI; yI]; IC = OC-OI; IC = IC/norm(IC); xN = IC(1); yN = IC(2); return;
@mirror/ray.m

c) La classe diaphragme

Le fichier diaphragm.m dans le dossier @diaphragm contient la définition de la classe diaphragm :


function obj=diaphragm(varargin)
% DIAPHRAGM is the constructor function for diaphragm object.
% The diaphragm class is primarily meant to be used to store optical characteristics of a diaphragm.
% obj=DIAPHRAGM(varargin)
% Called with no arguments, the constructor returns a default diaphragm object.
% Called with one argument that is a diaphragm object, the object is simply returned.
% Called with five arguments (name,Z,R,D), the constructor returns a new diaphragm object.
%  » name = name of the dioptre in the optical system
%  » Z = position of vertex on optical axis [m]
%  » R = radius of curvature [m]
%  » D = diameter of aperture [m]
%  « obj = diaphragm object
%
% © Eric Anterrieu (2004)
%
% See also: DIAPHRAGM/DISPLAY, DIAPHRAGM/GET.
%
switch (nargin)
case 0
  obj.name = '';
  obj.Z = 0;
  obj.R = 0;
  obj.D = 0;
  obj = class(obj,'diaphragm');
case 1
  if (isa(varargin{1},'diaphragm'))
    obj = varargin{1};
  else
    error('Wrong argument type.');
   end
case 4
  if (varargin{4} < 0), error('Diameter of diaphragm should be positive.');  end
  obj.name = varargin{1};
  obj.Z    = varargin{2};
  obj.R    = varargin{3};
  obj.D    = varargin{4};
  obj = class(obj,'diaphragm');
otherwise
  error('Wrong number of input arguments.');
 end
%
return;
@diaphragm/diaphragm.m

Dans le dossier @diaphragm, on trouve la sur-définition de fonctions pour afficher/modifier les valeurs des propriétés d’un objet de la classe diaphragm...


function display(obj)
% DISPLAY(obj) displays the optical characteristics of a diaphragm.
%  » obj = diaphragm object
%
% © Eric Anterrieu (2004)
%
% See also: DIAPHRAGM/DIAPHRAGM, DIAPHRAGM/GET, DIAPHRAGM/SET.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
disp(sprintf('             name of the diaphragm: %s',  obj.name));
disp(sprintf('position of vertex on optical axis: %g m',obj.Z));
disp(sprintf('               radius of curvature: %g m',obj.R));
disp(sprintf('          diameter of the aperture: %g',  obj.D));
%
return;
@diaphragm/display.m


function obj=set(obj,varargin)
% obj=SET(obj,property,value,...) sets diaphragm properties to specified values.
%  » obj = diaphragm object
%  » property = diaphragm property
%  » val = value of the property
% Valid properties are: 'Name', 'Position', 'Curvature', 'Aperture'
%
% © Eric Anterrieu (2004)
%
% See also: DIAPHRAGM/DIAPHRAGM, DIAPHRAGM/GET, DIAPHRAGM/DISPLAY.
%
property_argin = varargin;
while length(property_argin) >= 2
  property = property_argin{1};
  val = property_argin{2};
  property_argin = property_argin(3:end);
  switch property
  case 'Name'
    obj.name = val;
  case 'Position'
    obj.Z = val;
  case 'Curvature'
    obj.R = val;
  case 'Aperture'
    obj.D = val;
  otherwise
    error([property,' is not a valid diaphragm property.']);
  end
end
%
return;
@diaphragm/set.m


function val=get(obj,property)
% val=GET(obj,property) returns value of a diaphragm property.
%  » obj = diaphragm object
%  » property = diaphragm property
%  « val = value of the property
% Valid properties are: 'Name', 'Position', 'Curvature', 'Aperture'
%
% © Eric Anterrieu (2004)
%
% See also: DIAPHRAGM/DIAPHRAGM, DIAPHRAGM/SET, DIAPHRAGM/DISPLAY.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
ifisa(obj,'diaphragm')), error('Wrong argument type.');  end
ifisa(property,'char')), error('Wrong argument type.');  end
%
switch (property)
case 'Name'
  val = obj.name;
case 'Position'
  val = obj.Z;
case 'Curvature'
  val = obj.R;
case 'Aperture'
  val = obj.D;
otherwise
  error([property,' is not a valid diaphragm property.']);
end
%
return;
@diaphragm/get.m

Dans le dossier @diaphragm, on trouve aussi des fonctions internes à la classe diaphragm...


function mat=matrix(obj)
% mat=MATRIX(obj) returns the refraction matrix of a diaphragm.
%  » obj = diaphragm object
%  « mat = refraction matrix
%
% © Eric Anterrieu (2004)
%
% See also: DIAPHRAGM/DIAPHRAGM, DIAPHRAGM/GET.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
mat = eye(2);
%
return;
@diaphragm/matrix.m


function mat=propagation(obj1,obj2)
% mat=PROPAGATION(obj1,obj2) returns the propagation matrix from a diaphragm to another optical subsystem.
%  » obj1 = diaphragm object
%  » obj2 = optical subsystem
%  « mat = propagation matrix
%
% © Eric Anterrieu (2004)
%
% See also: DIAPHRAGM/DIAPHRAGM, DIAPHRAGM/MATRIX, DIAPHRAGM/GET.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
mat = eye(2);
if (isa(obj2,'diaphragm') == 1)
  error('Only one diaphragm between two dioptres/mirrors.');
elseif (isa(obj2,'dioptre') == 1)
  mat(1,2)=abs(get(obj2,'Position') - obj1.Z)/get(obj2,'FontIndex');
elseif (isa(obj2,'mirror') == 1)
  mat(1,2)=abs(get(obj2,'Position') - obj1.Z)/get(obj2,'FrontIndex');
end
%
return;
@diaphragm/propagation.m

Dans le dossier @diaphragm, on trouve aussi des fonctions destinées au tracé de rayons...


function h=plot(obj,clipping,varargin)
% h=PLOT(obj,clipping,varargin) draws a diaphragm and returns a handle to the line object.
%  » obj = diaphragm object
%  » clipping = optical subsystem
%  » varargin = line attributs
%  « h = handle to line object
%
% © Eric Anterrieu (2004)
%
% See also: PLOTSYSOPT, DIAPHRAGM/DIAPHRAGM.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
y=[-clipping -obj.D/2 NaN obj.D/2 clipping];
x=ones(size(y))*obj.Z;
if (nargout == 0),   plot(x,y,varargin{:}); end
if (nargout == 1), h=plot(x,y,varargin{:}); end
%
return;
@diaphragm/plot.m


function Rray=ray(obj,Iray)
% Rray=RAY(obj,Iray) returns the point hitted by an incidence ray on a diaphragm
% as well as the direction of the transmitted ray.
%  » obj = diaphragm object
%  » Iray(3) = incident ray
%  « Rray(3) = transmitted ray
% First and second components of the ray are the co-ordinates of a point on the ray,
% while the third one is the angle made by the ray with the optical axis.
%
% © Eric Anterrieu (2004)
%
% See also: DIAPHRAGM/DIAPHRAGM.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
if (isnan(Iray(3))), Rray=[NaN;NaN;NaN]; return; end
%
Z = obj.Z;
D = obj.D;
% I
a = tan(Iray(3));
b = Iray(2) - a*Iray(1);
Rray(1) = Z;
Rray(2) = a*Rray(1) + b;
% U2
if (abs(Rray(2)) <= D/2)
  Rray(3) = Iray(3);
else
  Rray(3) = NaN;
end
%
return;
@diaphragm/ray.m

I.2.- Matrice de transfert

La fonction matrix.m retourne la matrice de transfert T(ES) d’un système optique entre les plans d’entrée Exy et de sortie Sxy :


function T=matrix(sys)
% T=MATRIX(sys) returns the transfert matrix through an optical system
%  » sys = optical system
%  « T = 2x2 transfert matrix
%
% © Eric Anterrieu (2004)
%
% See also: SYSOPT, CARDINAL, FOCAL, VERGENCE.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
for k=2:length(sys)
  kk(k-1)=isa(sys{k},'diaphragm');
end
kk=1+find(kk == 0);
%
T=matrix(sys{2});
for k=3:length(sys)
  T = propagation(sys{k-1},sys{k})*T;
  T = matrix(sys{k})*T;
end
%
return;
matrix.m

La fonction vergence.m retourne la vergence V d’un système optique :


function V=vergence(sys)
% V=VERGENCE(sys) returns the power of an optical system
%  » sys = optical system
%  « V = power
%
% © Eric Anterrieu (2004)
%
% See also: SYSOPT, CARDINAL, FOCAL, MATRIX.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
T = matrix(sys)
V = -T(2,1);
if (nargout == 0)
  fprintf(1,'V   = %+.3f\n',V);
end
%
return;
vergence.m

I.3.- Éléments cardinaux

La fonction focal.m retourne les distances focales objet fo et image fi :


function [fo,fi]=focal(sys)
% [fo,fi]=FOCAL(sys) returns the focal lengths of an optical system
%  » sys = optical system
%  » fo = object focal length
%  « fi = image focal length
%
% © Eric Anterrieu (2004)
%
% See also: SYSOPT, CARDINAL, VERGENCE, MATRIX.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
V = vergence(sys);
if (V == 0)
  fo = -Inf;
  fi =  Inf;
  warning('Optical system is a non-focal one.');
else
  No =  sys{1}(1);
  Ni =  sys{1}(2);
  fo = -No/V;
  fi =  Ni/V;
end
if (nargout == 0)
  fprintf(1,'fo  = %+.3f m   fi  = %+.3f m\n',fo,fi);
end
%
return;
focal.m

La fonction cardinal.m renvoie les positions, par rapport aux plans d’entrée Exy et de sortie Sxy, des plans focaux EFo et SFi, des plans principaux EHo et SHI, et des points nodaux ENo et SNi :


function [EFo,SFi,EHo,SHi,ENo,SNi]=cardinal(sys)
% [EFo,SFi,EHo,SHi,ENo,SNi]=CARDINAL(sys) returns the position of the cardinal elements
%  » sys = optical system
%  « EFo,SFi = focal planes
%  « EHo,SHi = unit planes
%  « ENo,SNi = nodal points
%
% © Eric Anterrieu (2004)
%
% See also: PLOTCARDINAL, SYSOPT, FOCAL, VERGENCE, MATRIX.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
%
V = vergence(sys);
if (V == 0), error('Optical system is a non-focal one.'); end
%
T = matrix(sys);
[fo,fi] = focal(sys);
No = sys{1}(1);
Ni = sys{1}(2);
%
EFo = fo*T(2,2);
SFi = fi*T(1,1);
%
EHo = fo*(T(2,2)-1);
SHi = fi*(T(1,1)-1);
%
ENo = fo*(T(2,2)-Ni/No);
SNi = fi*(T(1,1)-No/Ni);
%
if (nargout == 0)
  fprintf(1,'EFo = %+.3f m  SFi = %+.3f m\n',EFo,SFi);
  fprintf(1,'EHo = %+.3f m  SHi = %+.3f m\n',EHo,SHi);
  fprintf(1,'ENo = %+.3f m  SNi = %+.3f m\n',ENo,SNi);
end
%
return;
cardinal.m

I.4.- Formules de conjugaison

Il existe trois relations de conjugaison permettant de calculer la position d’une image connaissant la position de l’objet, et vice versa. Elles diffèrent selon l’origine des points utilisés :

  • les points d’entrée et de sortie sur l’axe optique, pour la relation homographique ;
  • les points principaux, pour la relation de Descartes ;
  • les foyers, pour la relation de Newton.

a) Relation homographique

La relation homographique est mise en oeuvre dans les fonctions homographic.m, et respectivement invhomographic.m, afin de calculer la position d’une image Ai (respectivement d’un objet Ao) connaissant la position d’un objet Ao (respectivement d’une image Ai) par rapport au plan d’entré Exy (respectivement plan de sortie Sxy) du système optique. Par ailleurs, ces deux programmes renvoient aussi les grandissements transversal Gt, angulaire Ga et longitudinal Gl.


function [SAi,Gt,Ga,Gl]=homographic(sys,EAo)
% [SAi,Gt,Ga,Gl]=HOMOGRAPHIC(sys,EAo) returns the position of an image wrt. the input
% plane for a given object whose position is given wrt. the output plane
%  » sys = optical system
%  » EAo = position of the object on the optical axis
%  « SAi = position of the image on the optical axis
%  « Gt = transversal magnification
%  « Ga = angular magnification
%  « Gl = longitudinal magnification
%
% © Eric Anterrieu (2004)
%
% See also: MATRIX, DESCARTES, NEWTON, INVHOMOGRAPHIC, INVDESCARTES, INVNEWTON.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
T = matrix(sys);
No = sys{1}(1);
Ni = sys{1}(2);
SAi = Ni*(T(1,1)*(EAo/No)-T(1,2))/(-T(2,1)*(EAo/No)+T(2,2));
%
TAoAi = [1 SAi/Ni; 0 1]*T*[1 -EAo/No; 0 1];
Gt = TAoAi(1,1);
Ga = TAoAi(2,2)*No/Ni;
Gl = Gt/Ga;
%
return;
homographic.m


function [EAo,Gt,Ga,Gl]=invhomographic(sys,SAi)
% [EAo,Gt,Ga,Gl]=INVHOMOGRAPHIC(sys,SAi) returns the position of an image wrt. the output
% plane for a given object whose position is given wrt. the input plane
%  » sys = optical system
%  » SAi = position of the image on the optical axis
%  « EAo = position of the object on the optical axis
%  « Gt = transversal magnification
%  « Ga = angular magnification
%  « Gl = longitudinal magnification
%
% © Eric Anterrieu (2004)
%
% See also: MATRIX, INVDESCARTES, INVNEWTON, HOMOGRAPHIC, DESCARTES, NEWTON.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
T = matrix(sys);
No = sys{1}(1);
Ni = sys{1}(2);
EAo = No*(T(2,2)*(SAi/Ni)+T(1,2))/(T(2,1)*(SAi/Ni)+T(1,1));
%
TAoAi = [1 SAi/Ni; 0 1]*T*[1 -EAo/No; 0 1];
Gt = TAoAi(1,1);
Ga = TAoAi(2,2)*No/Ni;
Gl = Gt/Ga;
%
return;
invhomographic.m

b) Formule de Descartes

La relation de conjugaison de Descartes est mise en oeuvre dans les fonctions descartes.m, et respectivement invdescartes.m, afin de calculer la position d’une image Ai (respectivement d’un objet Ao) connaissant la position d’un objet Ao (respectivement d’une image Ai) par rapport au plan principal objet Ho (respectivement plan principal image Hi) du système optique. Par ailleurs, ces deux programmes renvoient aussi les grandissements transversal Gt, angulaire Ga et longitudinal Gl.


function [HiAi,Gt,Ga,Gl]=descartes(sys,HoAo)
% [HiAi,Gt,Ga,Gl]=DESCARTES(sys,HoAo) returns the position of an image wrt. the image principal
% plane for a given object whose position is given wrt. the object principal plane
%  » sys = optical system
%  » HoAo = position of the object on the optical axis
%  « HiAi = position of the image on the optical axis
%  « Gt = transversal magnification
%  « Ga = angular magnification
%  « Gl = longitudinal magnification
%
% © Eric Anterrieu (2004)
%
% See also: VERGENCE, NEWTON, HOMOGRAPHIC, INVDESCARTES, INVNEWTON, INVHOMOGRAPHIC.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
V = vergence(sys);
if (V == 0)
  HiAi = NaN;
  Gt = NaN;
  Ga = NaN;
  Gl = NaN;
  warning('Optical system is a non-focal one.');
else
  No = sys{1}(1);
  Ni = sys{1}(2);
  HiAi = Ni/(V + No/HoAo);
  Gt = (No/Ni)*(HiAi/HoAo);
  Ga = HoAo/HiAi;
  Gl = Gt/Ga;
end
%
return;
descartes.m


function [HoAo,Gt,Ga,Gl]=invdescartes(sys,HiAi)
% [HoAo,Gt,Ga,Gl]=DESCARTES(sys,HiAi) returns the position of an object wrt. the object principal
% plane for a given image whose position is given wrt. the image principal plane
%  » sys = optical system
%  » HiAi = position of the image on the optical axis
%  « HoAo = position of the object on the optical axis
%  « Gt = transversal magnification
%  « Ga = angular magnification
%  « Gl = longitudinal magnification
%
% © Eric Anterrieu (2004)
%
% See also: VERGENCE, INVNEWTON, INVHOMOGRAPHIC, DESCARTES, NEWTON, HOMOGRAPHIC.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
V = vergence(sys);
if (V == 0)
  HoAo = NaN;
  Gt = NaN;
  Ga = NaN;
  Gl = NaN;
  warning('Optical system is a non-focal one.');
else
  No = sys{1}(1);
  Ni = sys{1}(2);
  HoAo = No/(Ni/HiAi - V);
  Gt = (No/Ni)*(HiAi/HoAo);
  Ga = HoAo/HiAi;
  Gl = Gt/Ga;
end
%
return;
invdescartes.m

c) Formule de Newton

La relation de conjugaison de Newton est mise en oeuvre dans les fonctions newton.m, et respectivement invnewton.m, afin de calculer la position d’une image Ai (respectivement d’un objet Ao) connaissant la position d’un objet Ao (respectivement d’une image Ai) par rapport au plan focal objet Fo (respectivement plan focal image Fi) du système optique. Par ailleurs, ces deux programmes renvoient aussi les grandissements transversal Gt, angulaire Ga et longitudinal Gl.


function [FiAi,Gt,Ga,Gl]=newton(sys,FoAo)
% [FiAi,Gt,Ga,Gl]=NEWTON(sys,FoAo) returns the position of an image wrt. the image focal
% plane for a given object whose position is given wrt. the object focal plane
%  » sys = optical system
%  » FoAo = position of the object on the optical axis
%  « FiAi = position of the image on the optical axis
%  « Gt = transversal magnification
%  « Ga = angular magnification
%  « Gl = longitudinal magnification
%
% © Eric Anterrieu (2004)
%
% See also: VERGENCE, DESCARTES, HOMOGRAPHIC, INVNEWTON, INVDESCARTES, INVHOMOGRAPHIC.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
V = vergence(sys);
if (V == 0)
  FiAi = NaN;
  Gt = NaN;
  Ga = NaN;
  Gl = NaN;
  warning('Optical system is a non-focal one.');
else
  No = sys{1}(1);
  Ni = sys{1}(2);
  FiAi = -(No*Ni)/V^2/FoAo;
  Gt = -(V*FiAi)/Ni;
  Ga = -No/(V*FiAi);
  Gl = Gt/Ga;
end
%
return;
newton.m


function [FoAo,Gt,Ga,Gl]=invnewton(sys,FiAi)
% [FoAo,Gt,Ga,Gl]=INVNEWTON(sys,FiAi) returns the position of an object wrt. the object focal
% plane for a given image whose position is given wrt. the image focal plane
%  » sys = optical system
%  » FiAi = position of the image on the optical axis
%  « FoAo = position of the object on the optical axis
%  « Gt = transversal magnification
%  « Ga = angular magnification
%  « Gl = longitudinal magnification
%
% © Eric Anterrieu (2004)
%
% See also: VERGENCE, INVDESCARTES, INVHOMOGRAPHIC, NEWTON, DESCARTES, HOMOGRAPHIC.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
V = vergence(sys);
if (V == 0)
  FoAo = NaN;
  Gt = NaN;
  Ga = NaN;
  Gl = NaN;
  warning('Optical system is a non-focal one.');
else
  No = sys{1}(1);
  Ni = sys{1}(2);
  FoAo = -(No*Ni)/V^2/FiAi;
  Gt =  No/(V*FoAo);
  Ga =  (V*FoAo)/Ni;
  Gl = Gt/Ga;
end
%
return;
invnewton.m

I.5.- Construction géométrique dans l’approximation de Gauss

La fonction plotcardinal.m trace, sur un système optique représenté sous la forme schématique d’un domaine délimité par les plans d’entrée Exy et de sortie Sxy, la position des différents éléments cardinaux.


function plotcardinal(sys)
% PLOTCARDINAL(sys) draws 
%  » sys = optical system
%
% © Eric Anterrieu (2004)
%
% See also: CARDINAL, SYSOPT, FOCAL, VERGENCE, MATRIX.
%
if (nargin ˜= 1), error('Wrong number of input arguments.'); end
% non-focal system?
V=vergence(sys);
if (V == 0), error('Optical system is a non-focal one.'); end
% cardinal elements
[EFo,SFi,EHo,SHi,ENo,SNi]=cardinal(sys);
% E and S
E=sys{1}(3); S=sys{1}(4);
hf=fill([E E S S],[-1 1 1 -1],[1 1 1]*0.9); set(hf,'EdgeColor',get(hf,'FaceColor'));
plot(E,0,'k+',S,0,'k+');
ht=text(E,-0.05,'E'); set(ht,'Color','k');
ht=text(S,-0.05,'S'); set(ht,'Color','k');
% Fo and Fi
plot(E+EFo,0,'r*',S+SFi,0,'r*');
ht=text(E+EFo,-0.05,'F_o'); set(ht,'Color','r');
ht=text(S+SFi,-0.05,'F_i'); set(ht,'Color','r');
% Ho and Hi
plot(E+[1 1]*EHo,[-1 1],'m-',S+[1 1]*SHi,[-1 1],'m-');
ht=text(E+EHo,-1.05,'H_o'); set(ht,'Color','m');
ht=text(S+SHi,-1.05,'H_i'); set(ht,'Color','m');
% No and Ni
plot(E+ENo,0,'go',S+SNi,0,'go');
ht=text(E+ENo,0.05,'N_o'); set(ht,'Color','g');
ht=text(S+SNi,0.05,'N_i'); set(ht,'Color','g');
%
zmin=min([E+[EFo EHo ENo], S+[SFi SHi SNi]]);
zmax=max([E+[EFo EHo ENo], S+[SFi SHi SNi]]);
axis([zmin-0.03*(zmax-zmin) zmax+0.03*(zmax-zmin) -1.25 1.25]);
%
return;
plotcardinal.m

La fonction plotgaussrays trace, pour un objet de hauteur AoBo, situé à une distance EAo de Exy et pour une image de hauteur AiBi, située à une distance SAi de Sxy, les trois rayons suivants dans l’approximation de Gauss :

  • le rayon issu de Bo, parallèle à l’axe optique, qui émerge en passant par le foyer image Fi en direction de Bi ;
  • le ayon issu de Bo, passant par le foyer objet Fo, qui émerge parallèlement à l’axe optique en passant par Bi ;
  • le rayon issu de Bo, passant par le point nodal objet No, qui émerge vers Bi parallèlement à sa direction incidente en provenant du point nodal image Ni.

function plotgaussrays(sys,EAo,AoBo)
% PLOTGAUSSRAYS(sys,EAo,AoBo) draws 
%  » sys = optical system
%  » EAo = position of the object on the optical axis
%  « AoBo = size of the object
%
% © Eric Anterrieu (2004)
%
% See also: CARDINAL, SYSOPT, FOCAL, VERGENCE, MATRIX.
%
if (nargin ˜= 3), error('Wrong number of input arguments.'); end
% non-focal system?
V=vergence(sys);
if (V == 0), error('Optical system is a non-focal one.'); end
% cardinal elements
[EFo,SFi,EHo,SHi,ENo,SNi]=cardinal(sys);
plotcardinal(sys);
% image vs. object
[SAi,G] = homographic(sys,EAo);
AiBi = G*AoBo;
% AoBo
E=sys{1}(3);
plot(E+[1 1]*EAo,[0 AoBo],'k');
if (AoBo < 0)
  plot(E+EAo,AoBo,'kv');
  ht=text(E+EAo,0.05,'A_o'); set(ht,'Color','k');
  ht=text(E+EAo,AoBo-0.05,'B_o'); set(ht,'Color','k');
else
  plot(E+EAo,AoBo,'k^');
  ht=text(E+EAo,-0.05,'A_o'); set(ht,'Color','k');
  ht=text(E+EAo,AoBo+0.05,'B_o'); set(ht,'Color','k');
end
% AiBi
S=sys{1}(4);
plot(S+[1 1]*SAi,[0 AiBi],'k');
if (AiBi < 0)
  plot(S+SAi,AiBi,'kv');
  ht=text(S+SAi,0.05,'A_i'); set(ht,'Color','k');
  ht=text(S+SAi,AiBi-0.05,'B_i'); set(ht,'Color','k');
else
  plot(S+SAi,AiBi,'k^');
  ht=text(S+SAi,-0.05,'A_i'); set(ht,'Color','k');
  ht=text(S+SAi,AiBi+0.05,'B_i'); set(ht,'Color','k');
end
% rays from Bo to Bi
plotBoFoBi(sys,EAo,AoBo,SAi,AiBi);
plotBoFiBi(sys,EAo,AoBo,SAi,AiBi);
plotBoNoNiBi(sys,EAo,AoBo,SAi,AiBi);
%
zmin=min([E+[EAo EFo EHo ENo], S+[SAi SFi SHi SNi]]);
zmax=max([E+[EAo EFo EHo ENo], S+[SAi SFi SHi SNi]]);
axis([zmin-0.03*(zmax-zmin) zmax+0.03*(zmax-zmin) -1.25 1.25]);
%
return;

function plotBoFoBi(sys,EAo,AoBo,SAi,AiBi) % [EFo,SFi,EHo,SHi,ENo,SNi]=cardinal(sys); if (EAo == EFo), return; end % E=sys{1}(3); S=sys{1}(4); Co=min([E+[EAo EFo EHo ENo], S+[SAi SFi SHi SNi]]); Ci=max([E+[EAo EFo EHo ENo], S+[SAi SFi SHi SNi]]); % ao = -AoBo/(EFo-EAo); bo = -ao*(E+EFo); ai = 0.0; bi = ao*(E+EHo) + bo; % xE = E; yE = ao*E + bo; xCo = Co; yCo = ao*Co + bo; xAo = E+EAo; yAo = ao*(E+EAo) + bo; xHo = E+EHo; yHo = ao*(E+EHo) + bo; xFo = E+EFo; yFo = ao*(E+EFo) + bo; xS = S; yS = ai*S + bi; xCi = Ci; yCi = ai*Ci + bi; xAi = S+SAi; yAi = ai*(S+SAi) + bi; xHi = S+SHi; yHi = yHo; % if (EAo <= 0.0) % real AoBo plot([xFo xE xHo],[yFo yE yHo],'r:'); plot([xCo xE],[yCo yE],'r-'); else % virtual AoBo plot([xFo xE xHo],[yFo yE yHo],'r:'); plot([xAo xE],[yAo yE],'r:'); plot([xCo xE],[yCo yE],'r-'); end plot([xHo xHi],[yHo yHi],'r:'); if (SAi >= 0.0) % real AiBi plot([xS xHi],[yS yHi],'r:'); plot([xS xCi],[yS yCi],'r-'); else % virtual AiBi plot([xS xAi],[yS yAi],'r:'); plot([xS xHi],[yS yHi],'r:'); plot([xS xCi],[yS yCi],'r-'); end % return;
function plotBoFiBi(sys,EAo,AoBo,SAi,AiBi) % [EFo,SFi,EHo,SHi,ENo,SNi]=cardinal(sys); if (SFi == SHi), return; end % E=sys{1}(3); S=sys{1}(4); Co=min([E+[EAo EFo EHo ENo], S+[SAi SFi SHi SNi]]); Ci=max([E+[EAo EFo EHo ENo], S+[SAi SFi SHi SNi]]); % ao = 0.0; bo = AoBo; ai = -AoBo/(SFi-SHi); bi = -ai*(S+SFi); % xE = E; yE = ao*E + bo; xCo = Co; yCo = ao*Co + bo; xAo = E+EAo; yAo = ao*(E+EAo) + bo; xHo = E+EHo; yHo = ao*(E+EHo) + bo; xS = S; yS = ai*S + bi; xCi = Ci; yCi = ai*Ci + bi; xAi = S+SAi; yAi = ai*(S+SAi) + bi; xHi = S+SHi; yHi = yHo; xFi = S+SFi; yFi = ai*(S+SFi) + bi; % if (EAo <= 0.0) % real AoBo plot([xE xHo],[yE yHo],'r:'); plot([xCo xE],[yCo yE],'r-'); else % virtual AoBo plot([xE xHo],[yE yHo],'r:'); plot([xE xAo],[yE yAo],'r:'); plot([xCo xE],[yCo yE],'r-'); end plot([xHo xHi],[yHo yHi],'r:'); if (SAi >= 0.0) % real AiBi plot([xS xFi],[yS yFi],'r:'); plot([xS xHi],[yS yHi],'r:'); plot([xS xCi],[yS yCi],'r-'); else % virtual AiBi plot([xS xAi],[yS yAi],'r:'); plot([xS xFi],[yS yFi],'r:'); plot([xS xHi],[yS yHi],'r:'); plot([xS xCi],[yS yCi],'r-'); end % return;
function plotBoNoNiBi(sys,EAo,AoBo,SAi,AiBi) % [EFo,SFi,EHo,SHi,ENo,SNi]=cardinal(sys); if (EAo == ENo), return; end % E=sys{1}(3); S=sys{1}(4); Co=min([E+[EAo EFo EHo ENo], S+[SAi SFi SHi SNi]]); Ci=max([E+[EAo EFo EHo ENo], S+[SAi SFi SHi SNi]]); % ao = -AoBo/(ENo-EAo); bo = -ao*(E+ENo); ai = ao; bi = -ai*(S+SNi); % xE = E; yE = ao*E + bo; xCo = Co; yCo = ao*Co + bo; xAo = E+EAo; yAo = ao*(E+EAo) + bo; xNo = E+ENo; yNo = 0; xS = S; yS = ai*S + bi; xCi = Ci; yCi = ai*Ci + bi; xAi = S+SAi; yAi = ai*(S+SAi) + bi; xNi = S+SNi; yNi = 0; % if (EAo <= 0.0) % real AoBo plot([xE xNo],[yE yNo],'g:'); plot([xE xCo],[yE yCo],'g-'); else % virtual AoBo plot([xE xAo],[yE yAo],'g:'); plot([xE xNo],[yE yNo],'g:'); plot([xE xCo],[yE yCo],'g-'); end if (SAi >= 0.0) % real AiBi plot([xS xNi],[yS yNi],'g:'); plot([xS xCi],[yS yCi],'g-'); else % virtual AiBi plot([xS xAi],[yS yAi],'g:'); plot([xS xNi],[yS yNi],'g:'); plot([xS xCi],[yS yCi],'g-'); end % return;
plotgaussrays.m

I.6.- Tracé de rayons

La fonction plotsysopt.m trace les éléments constitutifs d’un système optique en invoquant les fonctions membres plot des classes dioptre, mirror et diaphragm :


function plotsysopt(sys,clipping)
% PLOTSYSOPT(sys,clipping) draws subsystems of an optical system
%  » sys = optical system
%  » clipping = clipping in the direction perpendicular to the optical axis
%
% © Eric Anterrieu (2004)
%
% See also: SYSOPT, DIOPTRE/PLOT, MIRROR/PLOT, DIAPHRAGM/PLOT.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
for k=2:length(sys)
  plot(sys{k},clipping,'b-','Linewidth',2);
end
%
return;
plotsysopt.m

La fonction rays.m retourne les coordonnées des points d’intersection successifs de rayons incidents avec les différents éléments constitutifs d’un système optique, au cours de leur propagation à travers ce dernier, en invoquant les fonctions membres ray des classes dioptre, mirror et diaphragm :


function [x,y,a]=rays(sys,x0,y0,a0,xe)
% [x,y,a]=RAYS(sys,x0,y0,a0,xe) returns the coordinates of the successive points hitted by rays
% that propagate through an optical system
%  » sys = optical system
%  » x0(m) = x coordinates of the incident rays
%  » y0(m) = y coordinates of the incident rays
%  » a0(m) = incidence angles of the incident rays
%  » ze = position of the screen on the optical axis
%  « x(m,n+2) = x coordinates of the successive points hitted by the rays
%  « y(m,n+2) = y coordinates of the successive points hitted by the rays
%  « a(m,n+2) = incidence angles of the successive points hitted by the rays
% m is the number of rays, n is the number of optical subsystems.
% If y0 is a vector and x0,a0 are scalars, the rays are parallel.
% If a0 is a vector and x0,y0 are scalars, the rays are divergent.
%
% © Eric Anterrieu (2004)
%
% See also: PLOTRAYS, DIOPTRE/RAY, MIRROR/RAY, DIAPHRAGM/RAY.
%
if (nargin ˜= 5), error('Wrong number of input arguments.'); end
%
if (length(y0) > 1)
  y(:,1) = y0(:);
  x(:,1) = x0*ones(size(y));
  a(:,1) = a0*ones(size(y))*pi/180;
elseif (length(a0) > 1)
  a(:,1) = a0(:)*pi/180;
  x(:,1) = x0*ones(size(a));
  y(:,1) = y0*ones(size(a));
else
  a(:,1) = a0*pi/180;
  x(:,1) = x0;
  y(:,1) = y0;
end
%
for k=1:size(a,1)
  for n=2:length(sys)
     newI = ray(sys{n},[x(k,n-1);y(k,n-1);a(k,n-1)]);
     x(k,n) = newI(1);
     y(k,n) = newI(2);
     a(k,n) = newI(3);
  end
  if isnan(a(k,n))
     x(k,n+1) = NaN;
     y(k,n+1) = NaN;
     a(k,n+1) = NaN;
  else
     x(k,n+1) = xe;
     y(k,n+1) = xe*tan(a(k,n))+ (y(k,n)-tan(a(k,n))*x(k,n));
     a(k,n+1) = a(k,n);
  end
end
%
return;
rays.m

La fonction plotrays.m trace le trajet des rayons lumineux calculés par la fonction rays.m :


function plotrays(x,y)
% PLOTRAYS(x,y) draws rays that propagate through an optical system
%  » x(:,:) = x coordinates of the successive points hitted by the rays
%  » y(:,:) = y coordinates of the successive points hitted by the rays
%
% © Eric Anterrieu (2004)
%
% See also: RAYS, DIOPTRE/RAY, MIRROR/RAY, DIAPHRAGM/RAY.
%
if (nargin ˜= 2), error('Wrong number of input arguments.'); end
%
for k=2:length(sys)
  plot(x(k,:),y(k,:),'r-','Linewidth',1);
end
%
return;
plotrays.m

II.- Systèmes optiques centrés

Avec les éléments précédents, on étudie dans cette partie quelques systèmes optiques en associant des objets des classes dioptre, mirror et diaphragm à l’aide de la fonction sysopt.m :

II.1.- Système dioptrique

Le script suivant étudie une lentille boule d’indice n=1,5 et de rayon r=1,5 cm travaillant dans l’air :


E = dioptre('E',0.0E-02, 1.5E-02,1.0,1.5);
S = dioptre('S',1.5E-02, inf    ,1.5,1.0);
soc = sysopt(E,S);
T = matrix(soc)
V = vergence(soc)
[fo,fi] = focal(soc)
[EFo,SFi,EHo,SHi,ENo,SNi] = cardinal(soc)
figure; clf; hold('on'); plotcardinal(soc);
%
E = dioptre('E',1.5E-02, inf    ,1.0,1.5);
S = dioptre('S',3.0E-02,-1.5E-02,1.5,1.0);
soc = sysopt(E,S);
T = matrix(soc)
V = vergence(soc)
[fo,fi] = focal(soc)
[EFo,SFi,EHo,SHi,ENo,SNi] = cardinal(soc)
figure; clf; hold('on'); plotcardinal(soc);
%
E = dioptre('E',0.0E-02, 1.5E-02,1.0,1.5);
S = dioptre('S',3.0E-02,-1.5E-02,1.5,1.0);
soc = sysopt(E,S);
T = matrix(soc)
V = vergence(soc)
[fo,fi] = focal(soc)
[EFo,SFi,EHo,SHi,ENo,SNi] = cardinal(soc)
figure; clf; hold('on'); plotcardinal(soc);
EAo=0.O; AoBo=0.2;
[SAi,Gt,Ga,Gl] = homographic(soc,EAo)
figure; clf; hold('on'); plotgaussrays(soc,EAo,AoBo);
[HiAi,Gt,Ga,Gl] = descartes(soc,EAo-EHo)
SAi=SHi+HiAi
[FiAi,Gt,Ga,Gl] = newton(soc,EAo-EFo)
SAi=SFi+FiAi
[x,y]=rays(soc,-2.0E-02,linspace(-1.0E-02,1.0E-02,30),0,5.0E-02);
figure; clf; hold('on'); plotsysopt(soc,1.5E-02); plotrays(x,y);
%
return;
scriptII1.m

II.2.- Système catadioptrique

Le script suivant étudie une lentille boule d’indice n=1,5 et de rayon r=1,5 cm dont l’une des faces est argentée :


E = dioptre('E',0.0E-02, 1.5E-02,1.0,1.5);
M =  mirror('M',3.0E-02,-1.5E-02,1.5);
S = dioptre('S',0.0E-02,-1.5E-02,1.5,1.0);
soc = sysopt(E,M,S);
T = matrix(soc)
TEM = propagation(E,M)*matrix(E)
TME = matrix(E)*propagation(M,E)
T = TME*matrix(M)*TEM
V = vergence(soc)
[fo,fi] = focal(soc)
[EFo,SFi,EHo,SHi] = cardinal(soc)
figure; clf; hold('on'); plotcardinal(soc);
%
R = -2/V
EC=invhomographic(sysopt(E),1.5E-02)
ES=invhomographic(sysopt(E),3.0E-02)
R = EC-ES
%
[x,y]=rays(soc,-2.0E-02,linspace(-1.0E-02,1.0E-02,30),0,EFo);
figure; clf; hold('on'); plotsysopt(soc,1.5E-02); plotrays(x,y);
%
return;
scriptII2.m

III.- L’oeil

Les fonctions legrand.m et gullstrand.m renvoient, tout comme la fonction sysopt.m, un vecteur d’objets des classes dioptre et diaphragm qui contiennent les caractéristiques optiques d’un oeil humain standard dans une attitude accommodée ou relaxée :


function eye=legrand(type,subtype)
% eye=LEGRAND(type,subtype) returns the optical elements of the LeGrand eye
%  » type = 'relaxed' or 'accomodated'
%  » subtype = 'all' or 'cornea' or 'lens'
%  « eye = LeGrand relaxed/accommodated eye
%
% © Eric Anterrieu (2004)
%
% See also: SYSOPT, GULLSTRAND, DIOPTRE, DIAPHRAGM.
%
if (strncmp(lower(type),'accomodated',1) == 1)
  cornea1 =   dioptre('cornea',0.00E-03, 7.80E-03,1.0000,1.3771);
  cornea2 =   dioptre('cornea',0.55E-03, 6.50E-03,1.3771,1.3374);
  iris    = diaphragm('iris'  ,3.60E-03, 6.00E-03,6.00E-03);
  lens1   =   dioptre('lens'  ,3.20E-03, 6.00E-03,1.3374,1.4270);
  lens2   =   dioptre('lens'  ,7.70E-03,-5.50E-03,1.4270,1.3360);
else
  cornea1 =   dioptre('cornea',0.00E-03, 7.80E-03,1.0000,1.3771);
  cornea2 =   dioptre('cornea',0.55E-03, 6.50E-03,1.3771,1.3374);
  iris    = diaphragm('iris'  ,3.60E-03,10.20E-03,6.00E-03);
  lens1   =   dioptre('lens'  ,3.60E-03,10.20E-03,1.3374,1.4200);
  lens2   =   dioptre('lens'  ,7.60E-03,-6.00E-03,1.4200,1.3360);
end
eye = sysopt(cornea1,cornea2,iris,lens1,lens2);
if (nargin == 2)
  if (strncmp(lower(subtype),'cornea',1) == 1)
    eye = sysopt(cornea1,cornea2);
  end
  if (strncmp(lower(subtype),'lens',1) == 1)
    eye = sysopt(lens1,lens2);
  end
end
%
return;
legrand.m


function eye=gullstrand(type,subtype)
% eye=GULLSTRAND(type,subtype) returns the optical elements of the Gullstrand eye
%  » type = 'relaxed' or 'accomodated'
%  » subtype = 'all' or 'cornea' or 'lens'
%  « eye = Gullstrand relaxed/accommodated eye
%
% © Eric Anterrieu (2004)
%
% See also: SYSOPT, LEGRAND, DIOPTRE, DIAPHRAGM.
%
if (strncmp(lower(type),'accomodated',1) == 1)
  cornea1 =   dioptre('cornea',0.000E-03, 7.760E-03,1.000,1.376);
  cornea2 =   dioptre('cornea',0.500E-03, 6.800E-03,1.376,1.336);
  iris    = diaphragm('iris'  ,3.600E-03, 5.333E-03,8.000E-03);
  lens1   =   dioptre('lens'  ,3.600E-03, 5.333E-03,1.336,1.386);
  lens2   =   dioptre('lens'  ,4.146E-03, 2.655E-03,1.386,1.406);
  lens3   =   dioptre('lens'  ,6.565E-03,-2.655E-03,1.406,1.386);
  lens4   =   dioptre('lens'  ,7.200E-03,-5.333E-03,1.386,1.336);
else
  cornea1 =   dioptre('cornea',0.000E-03, 7.760E-03,1.000,1.376);
  cornea2 =   dioptre('cornea',0.500E-03, 6.800E-03,1.376,1.336);
  iris    = diaphragm('iris'  ,3.600E-03,10.000E-03,8.000E-03);
  lens1   =   dioptre('lens'  ,3.600E-03,10.000E-03,1.336,1.386);
  lens2   =   dioptre('lens'  ,4.146E-03, 7.911E-03,1.386,1.406);
  lens3   =   dioptre('lens'  ,6.565E-03,-5.760E-03,1.406,1.386);
  lens4   =   dioptre('lens'  ,7.200E-03,-6.000E-03,1.386,1.336);
end
eye = sysopt(cornea1,cornea2,iris,lens1,lens2,lens3,lens4);
if (nargin == 2)
  if (strncmp(lower(subtype),'cornea',1) == 1)
    eye = sysopt(cornea1,cornea2);
  end
  if (strncmp(lower(subtype),'lens',1) == 1)
    eye = sysopt(lens1,lens2,lens3,lens4);
  end
end
%
return;
gullstrand.m

III.1.- Approximation de Gauss

Le script suivant effectue quelques calculs dans l’approximation de Gauss puis détermine la position des pupilles d’entrée Pe et de sortie Ps à l’aide de la relation homographique :


eye = legrand('relaxed');
T = matrix(eye)
[fo,fi] = focal(eye)
[EFo,SFi,EHo,SHi,ENo,SNi] = cardinal(eye)
V = vergence(eye)
cornea = legrand('relaxed','cornea');  V1 = vergence(cornea)
lens   = legrand('relaxed','lens');    V2 = vergence(lens)
e = 1.3374*(V1+V2-V)/(V1*V2)
%
[EPe,Gt] = invhomographic(cornea,3.05E-03)
[SPs,Gt] =    homographic(lens,  0.00E-03)
%
return;
scriptIII1.m

III.2.- Tracé de rayons et aberrations

Le script suivant trace la marche des rayons d’un pinceau de rayons incidents parallèles à l’axe optique, depuis un plan de front situé en amont de la cornée jusqu’à un plan de front situé au voisinage de la rétine. A partir des coordonnées des rayons qui émergent du cristallin et de leur inclinaison sur l’axe optique on représente les variations de la distance z à laquelle les rayons coupent l’axe optique en fonction du rayon r du pinceau incident (on observe qu’ils convergent d’autant plus qu’ils sont éloignés de l’axe optique), ainsi que les variations du diamètre d de la tache formée sur la rétine en fonction de la distance de convergence z (on observe que la tache de moindre diffusion se trouve dans dans un plan de front situé en amont du plan focal image donné pa l’approximation de Gauss) :


[x,y,alpha]=rays(eye,-2.0E-03,linspace(-3.0E-03,3.0E-03,20),0,26.0E-03);
figure; clf; hold('on');
plotsysopt(eye,5.0E-3);
plotrays(x,y);
axis([-2 26 -6 6]*1.0E-03);
xlabel('z  [mm]');
ylabel('x  [mm]');
%
a = tan(alpha(:,end));
b = y(:,end) - a.*x(:,end);
z = -b./a;
r = abs(y(:,1));
figure; clf;
plot(z*1.0E+03,r*1.0E+03,'r-');
axis([21 26 0 3]);
xlabel('z  [mm]');
ylabel('r  [mm]');
%
z = (21:0.01:26)*1E-03;
for k=1:length(z), d(k) = 2*max(a*z(k) + b); end
figure; clf;
plot(z*1.0E+03,d*1.0E+03,'r-');
axis([21 26 0 1.2]);
xlabel('z  [mm]');
ylabel('d  [mm]');
%
return;
scriptIII2.m

III.3.- Lazik

Le script suivant effectue la simulation d’une phtoablation cornéenne pour corriger un défaut de myopie. On représente notamment les variations du rayon de courbure R du dioptre antérieur de la cornée en fonction de l’épaisseur e de tissu cornéen brûlé sur l’axe optique, puis les variations du diamètre d de la tache de moindre diffusion dans le plan de la rétine toujours en fonction de e :


R = 7.25E-03;
cornea1 =   dioptre('cornea',0.00E-03,    R    ,1.0000,1.3771);
cornea2 =   dioptre('cornea',0.41E-03, 5.25E-03,1.3771,1.3374);
iris    = diaphragm('iris'  ,3.45E-03,10.33E-03,6.00E-03);
lens1   =   dioptre('lens'  ,3.45E-03,10.33E-03,1.3374,1.4200);
lens2   =   dioptre('lens'  ,7.60E-03,-6.17E-03,1.4200,1.3360);
eye = sysopt(cornea1,cornea2,iris,lens1,lens2);
%
D = 6.0E-03;
h = R - sqrt(R^2 - (D/2)^2);
e = (0:150)*1.0E-06;
R = 0.5*((h-e).^2 + (D/2)^2)./(h - e);
figure; clf;
plot(e*1.0E+06,R*1.0E+03,'r-');
axis([0 150 7 10]);
xlabel('\epsilon  [µm]');
ylabel('R  [mm]');
%
for k=1:length(e)
  eye{2} = set(eye{2},'Position',e(k),'Curvature',R(k));
  [x,y,alpha] = rays(eye,-2.0E-03,linspace(-D/2,D/2,50),0,25.0E-03);
  a = tan(alpha(:,end));
  b = y(:,end) - a.*x(:,end);
  d(k) = 2*max(a*24.8E-03 + b);
end
figure; clf;
plot(e*1.0E+06,d*1.0E+03,'r-');
axis([0 150 0 1.2]);
xlabel('\epsilon  [µm]');
ylabel('d  [mm]');
[dmin,k]=min(d);
[e(k)*1.0E+06, d(k)*1.0E+03, R(k)*1.0E+03]
%
return;
scriptIII3.m

IV.- Propagation dans les milieux non homogènes

On se propose ici de résoudre l’équation différentielle vectorielle décrivant la trajectoire d’un rayon lumineux dans un milieu non homogène.

IV.1.- Milieu bidimensionnel

Le fichier ode2D.m définit le problème que l’on cherche à résoudre :


function df=ode2D(l,f,ro)
% f = [x, y, px, py]
x = f(1); px = f(3);
y = f(2); py = f(4);
% EDO: dpx/dl, dpy/dl
dpxdl = 0.5*(-2*x*ro^2/(x^2+y^2)^2);
dpydl = 0.5*(-2*y*ro^2/(x^2+y^2)^2);
% df = [px, py, dpx/dl, dpy/dl]
df = [px; py; dpxdl; dpydl];
%
return;
ode2D.m

Le script suivant précise d’abord les conditions initiales, résout ensuite l’équation différentielle vectorielle du premier ordre décrite dans le fichier oed2D.m en invoquant le solveur ode45 qui utilise une méthode Runge-Kutta d’ordre 5, puis représente enfin graphiquement la trajectoire du rayon solution :


% domain
ro = 1;
R  = 4*ro;
Ni = sqrt(1 + (ro/R)^2);
% x and y at starting point
h  = 2*ro;
Xi = sqrt(R^2 - h^2);
Yi = h;
% dx/dl and dy/dl at starting point
alpha = pi;
dXidl = Ni*cos(alpha);
dYidl = Ni*sin(alpha);
% solve ODE
options = odeset('Stats','on','RelTol',1E-12,'AbsTol',1E-12);
[l,f] = ode45('ode2D',[0:0.1:10],[Xi Yi dXidl dYidl],options,ro);
X = f(:,1);  Y = f(:,2);  dXdl = f(:,3);  dYdl = f(:,4);
% last point
k = max(find((X.^2+Y.^2) < R^2));
le    = interp1([X(k)^2+Y(k)^2,X(k+1)^2+Y(k+1)^2],[l(k),l(k+1)],R^2,'linear');
Xe    = interp1([l(k),l(k+1)],[X(k),X(k+1)],le,'linear');
Ye    = interp1([l(k),l(k+1)],[Y(k),Y(k+1)],le,'linear');
dXedl = interp1([l(k),l(k+1)],[dXdl(k),dXdl(k+1)],le,'linear');
dYedl = interp1([l(k),l(k+1)],[dYdl(k),dYdl(k+1)],le,'linear');
beta = atan2(dYedl,dXedl);
chi  = (pi+beta);
% trajectory
XA = Xi-R*cos(alpha);  XB = Xe+R*cos(beta);
YA = Yi-R*sin(alpha);  YB = Ye+R*sin(beta);
figure; clf; hold('on');
fill(4*cos((0:360)*pi/180),4*sin((0:360)*pi/180),[1 1 1]*0.9);
plot(  cos((0:360)*pi/180),  sin((0:360)*pi/180),'b-');
plot([XA; X(1:k); Xe; XB]/ro,[YA; Y(1:k); Ye; YB]/ro,'r-');
axis([-6 6 -6 6]);
xlabel('x/\rho_o');
ylabel('y/\rho_o');
%
return;
scriptIV1.m

IV.2.- Milieu tridimensionnel

Le fichier ode3D.m définit le problème que l’on cherche à résoudre :


function df=ode3D(l,f,a,Nc,D)
% f = [x, y, z, px, py, pz]
x = f(1); px = f(4);
y = f(2); py = f(5);
z = f(3); pz = f(6);
% EDO: dpx/dl, dpy/dl, dpz/dl
dpxdl = 0.5*(-4*Nc^2*D*x/a^2);
dpydl = 0.5*(-4*Nc^2*D*y/a^2);
dpzdl = 0;
% df = [px, py, pz, dpx/dl, dpy/dl, dpz/dl]
df = [px; py; pz; dpxdl; dpydl; dpzdl];
%
return;
ode3D.m

Le script suivant précise d’abord les conditions initiales, résout ensuite l’équation différentielle vectorielle du premier ordre décrite dans le fichier oed3D.m en invoquant le solveur ode45 qui utilise une méthode Runge-Kutta d’ordre 5, puis représente enfin graphiquement la trajectoire du rayon solution :


% domain
a = 1;
D = 0.01;
Nc = 1.46;
% x, y and z at starting point
ri = 0*a;
ti = 0*pi/180;
Xi = ri*cos(ti);
Yi = ri*sin(ti);
Zi = 0;
% N at starting point
Ni = Nc*sqrt(1 - 2*D*(ri/a)^2);
% dx/dl, dy/dl and dz/dl at starting point
alphaX = 2*pi/180;
alphaY = 0*pi/180;
dXidl = Ni*(sin(alphaX)*cos(ti) - cos(alphaX)*sin(alphaY)*sin(ti));
dYidl = Ni*(sin(alphaX)*sin(ti) + cos(alphaX)*sin(alphaY)*cos(ti));
dZidl = Ni*(cos(alphaX)*cos(alphaY));
% solve ODE
options = odeset('Stats','on','RelTol',1E-12,'AbsTol',1E-12);
[l,f] = ode45('ode3D',[0:0.5:100],[Xi Yi Zi dXidl dYidl dZidl],options,a,Nc,D);
X = f(:,1); Y = f(:,2); Z = f(:,3); dXdl = f(:,4); dYdl = f(:,5); dZdl = f(:,6);
% trajectory
figure; clf; hold('on');
plot3(Y,Z,X,'r-');
plot3(-0.8*ones(size(Y)),Z,X,'b-');
plot3(Y,160*ones(size(Z)),X,'b-');
plot3(Y,Z,-0.8*ones(size(X)),'b-');
axis([-0.8 0.8 0 160 -0.8 0.8]);
view(55,20); set(gca,'DataAspectRatio',[1 40 1]);
xlabel('y/a');
ylabel('z/a');
zlabel('x/a');
%
return;
scriptIV2.m

V.- Faisceaux gaussiens

On s’intéresse ici à la transformation d’un faisceau gaussien par un système optique. On se place dans le cas d’un faisceau monochromatique, issu d’un laser He-Ne, de longueur d’onde 488 nm et de waist 0,16 mm.

V.1.- Système dioptrique

Le script suivant étudie une lentille boule d’indice n=1,5 et de rayon r=1,5 cm travaillant dans l’air :


E = dioptre('E',0.0E-02, 1.5E-02,1.0,1.5);
S = dioptre('S',3.0E-02,-1.5E-02,1.5,1.0);
soc = sysopt(E,S);
T = matrix(soc)
%
EAo = -0.O9;
SAi = homographic(soc,EAo)
%
lambda = 488E-09;
wo = 0.16E-03;
Qo = EAo + j*(pi*wo^2/lambda);
Qi = homographic(soc,Qo);
SAi = real(Qi)
wi = sqrt(imag(Qi)/pi*lambda)
wi/wo
%
return;
scriptV1.m

V.2.- Système catadioptrique

Le script suivant étudie une lentille boule d’indice n=1,5 et de rayon r=1,5 cm dont l’une des faces est argentée :


E = dioptre('E',0.0E-02, 1.5E-02,1.0,1.5);
M =  mirror('M',3.0E-02,-1.5E-02,1.5);
S = dioptre('S',0.0E-02,-1.5E-02,1.5,1.0);
soc = sysopt(E,M,S);
T = matrix(soc)
%
EAo = -0.O9;
SAi = homographic(soc,EAo)
%
lambda = 488E-09;
wo = 0.16E-03;
Qo = EAo + j*(pi*wo^2/lambda);
Qi = homographic(soc,Qo);
SAi = real(Qi)
wi = sqrt(imag(Qi)/pi*lambda)
wi/wo
%
return;
scriptV2.m

VI.- Facteurs de visibilité en optique ondulatoire

On retrouve, par la simulation numérique, la distribution d’intensité, dans le plan image, à partir du facteur de visibilité, d’abord en éclairage incohérent puis en éclairage cohérent.

VI.1.- Facteur de visibilité en éclairage incohérent


L   = 597E-09;
I1  = 1.0;
mu  = 1.5;
DTo = 0.01E-03;
%
a = (0:1:200)*1.0E-03;
G = (exp(j*pi*DTo*a/L) + mu*exp(-j*pi*DTo*a/L))/(1+mu);
figure(1); clf; hold('on'); plot(a*1000,abs(G),'r-'); axis([0 200 0 1.02]);
figure(2); clf; hold('on'); plot(a*1000,angle(G),'g-'); axis([0 200 -1.05*pi 1.05*pi]);
%
a  = [5 18 30]*1.0E-03; 
Ti = 0:0.1*DTo:20*DTo;
for k=1:length(a)
  G  = (exp(j*pi*DTo*a(k)/L) + mu*exp(-j*pi*DTo*a(k)/L))/(1+mu);
  Ii = I1*(1+mu)*(1 + real(G*exp(j*2*pi*Ti*a(k)/L)));
  fprintf(1,'a=%.1f mm: max=%.5f min=%.5f V=%.5f mean=%.5f\n',
             a(k)*1000,max(Ii),min(Ii),(max(Ii)-min(Ii))/(max(Ii)+min(Ii)),mean(Ii));
  figure(1); plot(a(k)*1000,abs(G),'bo');
  figure(2); plot(a(k)*1000,angle(G),'bo');
  figure(k+2);
  subplot(2,1,1); imagesc(Ii,[0 5]); colormap('gray'); caxis([0 5]);
  subplot(2,1,2); plot(Ti*1000,Ii,'b-'); axis([0 0.2 0 5]);
  clear Ii G;
end
clear Ti;
%
D  = 50.0E-03;
To = -4*DTo:0.02*DTo:4*DTo;
Io = I1*( (sinc((To-0.5*DTo)*D/L)).^2 + mu*(sinc((To+0.5*DTo)*D/L)).^2 );
k1 = min(find(Io > 0.5*max(Io)));
k2 = max(find(Io > 0.5*max(Io)));
fwhm = To(k2)-To(k1);
fprintf(1,'D=%.1f mm: fwhm=%.2E\n',D*1000,fwhm);
figure(6); plot(To*1000,Io,'m-');
clear To;
%
a  = [10 15 20 25]*1.0E-03; 
Ti = 0:0.1*DTo:20*DTo;
for k=1:length(a)
  G  = (exp(j*pi*DTo*a(k)/L) + mu*exp(-j*pi*DTo*a(k)/L))/(1+mu);
  Ii = I1*(1+mu)*(1 + real(G*exp(j*2*pi*Ti*a(k)/L)));
  V(k)=(max(Ii)-min(Ii))/(max(Ii)+min(Ii));   
  fprintf(1,'a=%.1f mm: mu=%.5f V=%.5f\n',a(k)*1000,mean(Ii)-1,V(k));
end
%
sol = fmins('Finc',[1.53 0.0186E-03],[],'',a,V,L);
fprintf(1,'  µ = %.5f\n',sol(1));
fprintf(1,' Dx = %.2E rad\n',sol(2));
%
return;
scriptVI1.m

Le fichier Finc.m met en oeuvre l’équation non linéaire que l’on cherche à satisfaire, au sens des moindres carrés, dans la dernière partie du script précédent où l’on veut restaurer l’objet à partir du facteur de visibilité Vinc :


function f=Finc(x,a,V,L)
%
mu  = x(1);
DTo = x(2);
%
Vinc = sqrt(1 + mu^2 + 2*mu*cos(2*pi*DT0*a/L))/(1+mu);
f = sum((V-Vinc).^2);
%
return;
Finc.m

VI.2.- Facteur de visibilité en éclairage cohérent


L   = 597E-09;
I1  = 1.0;
mu  = 1.5;
DTo = 0.01E-03;
phi = 0.8;
%
a = (0:1:200)*1.0E-03;
G = (exp(j*pi*DTo*a/L) + mu*exp(-j*pi*DTo*a/L) + 2*sqrt(mu)*cos(phi));
G = G./(1 + mu + 2*sqrt(mu)*cos(phi)*cos(pi*DTo*a/L));
figure(1); clf; hold('on'); plot(a*1000,abs(G),'r-'); axis([0 200 0 1.02]);
figure(2); clf; hold('on'); plot(a*1000,angle(G),'g-'); axis([0 200 -1.05*pi 1.05*pi]);
%
a  = [5 33 44]*1.0E-03; 
Ti = 0:0.1*DTo:20*DTo;
for k=1:length(a)
  G = (exp(j*pi*DTo*a(k)/L) + mu*exp(-j*pi*DTo*a(k)/L) + 2*sqrt(mu)*cos(phi));
  G = G./(1 + mu + 2*sqrt(mu)*cos(phi)*cos(pi*DT0*a(k)/L));
  Ii = I1*(1+mu)*(1 + real(G*exp(j*2*pi*Ti*a(k)/L)));
  fprintf(1,'a=%.1f mm: max=%.5f min=%.5f V=%.5f mean=%.5f\n',
             a(k)*1000,max(Ii),min(Ii),(max(Ii)-min(Ii))/(max(Ii)+min(Ii)),mean(Ii));
  figure(1); plot(a(k)*1000,abs(G),'bo');
  figure(2); plot(a(k)*1000,angle(G),'bo');
  figure(k+2);
  subplot(2,1,1); imagesc(Ii,[0 5]); colormap('gray'); caxis([0 5]);
  subplot(2,1,2); plot(Ti*1000,Ii,'b-'); axis([0 0.2 0 5]);
  clear Ii G;
end
clear Ti;
%
D  = 50.0E-03;
To = -4*DTo:0.02*DTo:4*DTo;
Io = I1*( sinc((To-0.5*DTo)*D/L) + sqrt(mu)*sinc((To+0.5*DTo)*D/L) ).^2;
k1 = min(find(Io > 0.5*max(Io)));
k2 = max(find(Io > 0.5*max(Io)));
fwhm = To(k2)-To(k1);
fprintf(1,'D=%.1f mm: fwhm=%.2E\n',D*1000,fwhm);
figure(6); plot(To*1000,Io,'m-');
clear To;
%
a  = [10 14 18 22 26 30]*1.0E-03; 
Ti = 0:0.1*DTo:20*DTo;
for k=1:length(a)
  G = (exp(j*pi*DTo*a(k)/L) + mu*exp(-j*pi*DTo*a(k)/L) + 2*sqrt(mu)*cos(phi));
  G = G./(1 + mu + 2*sqrt(mu)*cos(phi)*cos(pi*DT0*a(k)/L));
  Ii = I1*(1+mu)*(1 + real(G*exp(j*2*pi*Ti*a(k)/L)));
  V(k)=(max(Ii)-min(Ii))/(max(Ii)+min(Ii));   
  fprintf(1,'a=%.1f mm: mu=%.5f V=%.5f\n',a(k)*1000,mean(Ii)-1,V(k));
end
%
sol = fmins('Fcoh',[1.52 0.0124E-03 pi/2],[],'',a,V,L);
fprintf(1,'  µ = %.5f\n',sol(1));
fprintf(1,' Dx = %.2E rad\n',sol(2));
fprintf(1,'phi = %.5f rad\n',sol(3));
%
return;
scriptVI2.m

Le fichier Fcoh.m met en oeuvre l’équation non linéaire que l’on cherche à satisfaire, au sens des moindres carrés, dans la dernière partie du script précédent où l’on veut restaurer l’objet à partir du facteur de visibilité Vcoh :


function f=Fcoh(x,a,V,L)
%
mu  = x(1);
DTo = x(2);
phi = x(3);
%
Vcoh = sqrt(1 + mu^2 + 2*mu*cos(2*pi*DT0*a/L) ...
              + 4*mu*(cos(phi))^2 + 4*sqrt(mu)*(1+mu)*cos(phi)*cos(pi*DT0*a/L)) ...
           ./(1 + mu + 2*sqrt(mu)*cos(phi)*cos(pi*DT0*a/L));
f = sum((V-Vcoh).^2);
%
return;
Fcoh.m

 

 

[Page Précédente] [Dans la même rubrique] [Sommaire]

 

 

Site du CNRS Logo OMP
Logo UPS
Logo INSU-CNRS

Annuaire

Rechercher

Sur ce site



Bibliotheque

Laboratoire Astrophysique de Toulouse - Tarbes (UMR5572)

CNRS (Midi-Pyrénées)

Univ. Paul Sabatier

Liens utiles

Accueil Imprimer Contact mail Plan du site Crédits