Fourth-Order Wavefront Aberrations

ford.m

function [W, L, epsilon0] = ford(pya, pyc, cv, th, rn, lambda, abbe)
%
% FORD
%
%  Fourth-order wavefront aberrations
%
%

deln = zeros(size(rn));
if (nargin>6)
    idx = find(abbe);
    deln(idx) = (rn(idx)-1)./abbe(idx);
end

m = length(th);
epsilon0 = abs(lambda/(rn(m)*pya(m,2)));
L = rn(m).*(pya(m,1).*pyc(m,2)-pya(m,2).*pyc(m,1));
for (i=2:m)
  Ia(i) = rn(i)*(pya(i,2)+cv(i)*pya(i,1));
  Ic(i) = rn(i)*(pyc(i,2)+cv(i)*pyc(i,1));
  Ja(i) = pya(i,1)*(pya(i,2)/rn(i) - pya(i-1,2)/rn(i-1));
  Jc(i) = pya(i,1)*(pyc(i,2)/rn(i) - pyc(i-1,2)/rn(i-1));
  Jp(i) = cv(i)*(1.0/rn(i)-1.0/rn(i-1));
  Jw(i) = deln(i)/rn(i) - deln(i-1)/rn(i-1);
end
W040 = -Ia.*Ia.*Ja/(8*lambda);
W131 = -Ia.*Ic.*Ja/(2*lambda);
W222 = -Ic.*Ic.*Ja/(2*lambda);
W220 = -L*L*Jp/(4*lambda) + 0.5*W222;
W311 = -Ic.*(Ic.*Jc-pyc(1:m,1)'.*Jp*L)/(2*lambda);
WchA = Ia.*pya(1:m,1)'.*Jw/(2*lambda);
WchL = Ic.*pya(1:m,1)'.*Jw/lambda;
W = [W040' W131' W222' W220' W311' WchA' WchL'];
W(m+1,:) = sum(W);

Example

Define the optical system.

>> cv = [ 0 1/34.53 -1/21.98 -1/214.63];
>> th = [ 0 9 2.5 0 ];
>> rn = [ 1 1.67 1.728 1 ];
>> abbe = [0 47.1 28.4 0];

Trace the axial ray, insert a thickness solve for the last surface to find the paraxial image plane, and repeat the raytrace.

>> pya = parax([6.25 0],cv,th,rn);
>> th(4)=-pya(4,1)/pya(4,2);
>> pya = parax([6.25 0],cv,th,rn);
>> pya

pya =

    6.2500         0
    6.2500   -0.0726
    5.5964   -0.0616
    5.4424   -0.1250
         0   -0.1250

Trace a temporary chief ray directed at the center of the first surface. Then modify it to make the stop at the last surface.

>> pyc = parax([0 tan(3.5*pi/180)],cv,th,rn);
>> pyc

pyc =

         0    0.0612
         0    0.0366
    0.3296    0.0359
    0.4194    0.0606
    3.0590    0.0606

>> k = -pyc(4,1)/pya(4,1);
>> k

k =

   -0.0771

>> pyc = pyc+k*pya;
>> pyc

pyc =

   -0.4816    0.0612
   -0.4816    0.0422
   -0.1016    0.0406
         0    0.0702
    3.0590    0.0702

Define the wavelength and calculate the aberrations.

>> lambda = 0.00055;
>> [W L e0] = ford(pya, pyc, cv, th, rn, lambda, abbe);
>> W

W =

         0         0         0         0         0         0         0
    2.0236    2.1114    0.5508    1.0471    0.5463    8.7601    4.5702
   -2.9688    1.6999   -0.2433   -0.1824    0.0522  -17.5615    5.0278
    2.4957   -4.6647    2.1796    1.2202   -1.1403   11.0327  -10.3103
    1.5506   -0.8533    2.4871    2.0849   -0.5418    2.2313   -0.7123
>> L

L =

    0.3823

>> e0

e0 =

    0.0044

Maintained by John Loomis, last updated 28 Nov 2001