Trigonometric Raytrace

Matlab code

function Y = trig_ray(Qin,Uin,cv,th,rn)
%
%  TRIG_RAY  trigonometric raytrace
%
m = length(th);

Q = Qin;
U = Uin;

for (n=2:m)
   I = asin(sin(U)+Q*cv(n));
   Iprime = asin(rn(n-1)*sin(I)/rn(n));
   Uprime = U + Iprime - I;
   G = Q./(cos(U)+cos(I));
   Qprime = G.*(cos(Uprime)+cos(Iprime));
   beta = Iprime - Uprime;
   Y = G.*(1+cos(beta));
   % uncomment following for intermediate results
   % [n Y U I Iprime Q Qprime G]
   Q = Qprime + th(n)*sin(Uprime);
   U = Uprime;
end

Example

First we define the optical system.

>> cv = [ 0 0.034654 -0.005408 0];
>> th = [ 0 1 49.42301 0];
>> rn = [ 1 1.5 1 1 ];

Then we can generate a single trigonmetric ray:

>> format short g
>> ye = trig_ray(6.25, 0, cv, th, rn);

ye = -0.10864

Meridional Ray Fan

We can use trig_ray to trace a fan of rays:

>> y = linspace(0,1,6)';
>> ye = trig_ray(y*6.25,0,cv,th,rn);
>> [y ye]
ans =
            0            0
          0.2  -0.00083418
          0.4   -0.0067075
          0.6     -0.02283
          0.8    -0.054764
            1     -0.10864

We can fit the wavefront polynomial as follows:

>> lambda = 0.00055;
>> ua = -0.125;
>> epsilon = abs(lambda/ua);

>> d = - ye/epsilon;
>> c = [ 2*y 4*y.^3 6*y.^5 ];
>> a = c\d

a =
    0.0022615
       5.9053
      0.17763

The rms deviation can be calculated as:

>> diff = d - c*a;
>> rms = sqrt(sum(diff.*diff))/length(diff)

rms =
   0.00018571

Maintained by John Loomis, last updated 24 May 1999