function [ pos, dir, q ] = raytrace(pos_in, dir_in, cv, th, rn)
%
% RAYTRACE
%
nsurf = length(th);
opl = 0.0;
pos = pos_in;
dir = dir_in;
fprintf('%4s%10s%10s%10s%8s%8s%8s%12s\n','n','X','Y','Z','kx','ky','kz','OPL');
fprintf('%4d%10.3f%10.3f%10.3f%8.4f%8.4f%8.4f\n',0,pos,dir);
for (n=2:nsurf-1)
[pos dir q(n)] = surf_trace(pos,dir,cv(n),th(n-1),rn(n),rn(n-1));
fprintf('%4d%10.3f%10.3f%10.3f%8.4f%8.4f%8.4f%12.8f\n',n-1,pos,dir,q(n));
opl = opl + rn(n-1)*rn(n-1)*q(n);
end
[pos dir qt] = surf_trace(pos,dir,cv(nsurf),th(nsurf-1),rn(nsurf),rn(nsurf-1));
fprintf('%4d%10.3f%10.3f%10.3f%8.4f%8.4f%8.4f\n',n,pos,dir);
Trace ray through single surface (used by raytrace)
function [pos, dir, q] = surf_trace(pos_in, dir_in, cv, th, rno, rni)
% SURF_TRACE trace ray through single spherical surface
% [POS DIR Q] = SURF_TRACE(POS_IN, DIR_IN, CV, TH, RNO, RNI)
% CV surface curvature
% TH axial distance to current surface
% RNO index of refraction after surface
% RNI index of refraction before surface
% POS_IN initial ray position at previous surface
% DIR_IN optical direction cosines of initial ray
%
% Output variables:
% POS ray position at current surface
% DIR optical direction cosines after refraction
% Q reduced distance to current surface
%
pos = pos_in;
pos(3) = pos(3) - th;
% intersect current surface
if (cv ~= 0)
A = cv*rni*rni;
B = dir_in(3) - cv*dot(pos,dir_in);
C = cv*dot(pos,pos)-2.0*pos(3);
D = B*B-A*C;
%fprintf('Quadratic (A = %g, B = %g, C = %g, D = %g)\n',A,B,C,D);
D = sqrt(D);
if (B<0)
D = -D;
end
q = C/(B+D);
else
q = -pos(3)/dir_in(3);
end
% calculate point of intersection
pos = pos + q*dir_in;
% calculate surface normal
norm = -cv*pos;
norm(3) = 1.0 + norm(3);
%fprintf('surface normal %g %g %g\n',norm);
% refract ray into surface
gcosi = dot(dir_in,norm);
root = gcosi*gcosi + (rno+rni)*(rno-rni);
gcosr = sqrt(root);
if (gcosi<0.0)
gcosr = -gcosr;
end
% calculate refracted ray direction
power = gcosr - gcosi;
dir = dir_in + power * norm;
% calculate optical path length
q = q * rni*rni;
First we define the optical system.
% define optical system cv = [ 0 1/34.53 -1/21.98 -1/214.63 0]; th = [ 0 9 2.5 43.567826 0]; rn = [ 1 1.670028 1.72825 1 1 ];
Then we define the starting ray
pos = [0 12.5 0]; dir = [0 0 1];
Call raytrace to do the work
[pos dir q] = raytrace(pos,dir,cv,th,rn);
This produces
n X Y Z kx ky kz OPL 0 0.000 12.500 0.000 0.0000 0.0000 1.0000 1 0.000 12.500 2.342 0.0000 -0.2527 1.6508 2.34194476 2 0.000 12.029 -3.584 0.0000 -0.2105 1.7154 5.19366516 3 0.000 11.319 -0.299 0.0000 -0.2500 0.9683 10.07334993 4 0.000 -0.005 0.000 0.0000 -0.2500 0.9683
Display the intersection of the ray with the image plane
fprintf('image plane: (y = %.6f, x = %.6f)',pos(2),pos(1));
The final result is
image plane: (y = -0.004580, x = 0.000000)
Add a non-refracting reference surface as the last surface in the optical system
cv = [ 0 1/34.53 -1/21.98 -1/214.63 1/43.567826 0]; th = [ 0 9 2.5 0 43.567826 0]; rn = [ 1 1.670028 1.72825 1 1 1 ];
Trace a ray through the edge of the aperture, as before. This time display the values of the optical pathlength array q.
pos = [0 12.5 0];
dir = [0 0 1];
[pos dir q] = raytrace(pos,dir,cv,th,rn);
fprintf('image plane: (y = %.6f, x = %.6f)\n',pos(2),pos(1));
disp(q);
The results are
n X Y Z kx ky kz OPL 0 0.000 12.500 0.000 0.0000 0.0000 1.0000 1 0.000 12.500 2.342 0.0000 -0.2527 1.6508 2.34194476 2 0.000 12.029 -3.584 0.0000 -0.2105 1.7154 5.19366516 3 0.000 11.319 -0.299 0.0000 -0.2500 0.9683 10.07334993 4 0.000 10.886 1.382 0.0000 -0.2500 0.9683 1.73561739 5 0.000 -0.005 0.000 0.0000 -0.2500 0.9683 image plane: (y = -0.004580, x = 0.000000) q = 0 2.3419 5.1937 10.0733 1.7356
Now trace a central ray. The optical pathlength array is named qc.
pos = [ 0 0 0]; dir = [ 0 0 1]; [pos dir qc] = raytrace(pos,dir,cv,th,rn); disp(qc);
The results are
n X Y Z kx ky kz OPL 0 0.000 0.000 0.000 0.0000 0.0000 1.0000 1 0.000 0.000 0.000 0.0000 0.0000 1.6700 0.00000000 2 0.000 0.000 -0.000 0.0000 0.0000 1.7283 15.03025200 3 0.000 0.000 0.000 0.0000 0.0000 1.0000 4.32062500 4 0.000 0.000 0.000 0.0000 0.0000 1.0000 0.00000000 5 0.000 0.000 0.000 0.0000 0.0000 1.0000 qc = 0 0 15.0303 4.3206 0
The final optical path is the sum of the differences in optical path length
opd = qc - q
opd = sum(opd);
fprintf('total opd %g\n',opd);
wvl = 0.00058756;
fprintf('total opd %g wave\n',opd/wvl);
The final results are
opd = 0 -2.3419 9.8366 -5.7527 -1.7356 total opd 0.00629975 total opd 10.7219 wave
Maintained by John Loomis, last updated 6 Feb 2003