general raytrace
Contents
Surface trace
cv = 1/34.53; % surface curvature rno=1.670028; % index of refraction after surface (out) rni=1; % index of refraction before surface (in) th=0; % distance from last surface pos_in = [0 3 -10]; % starting position of ray dir_in = [0,0.5,1]; % direction of ray (not normalized) dir_in = dir_in/norm(dir_in); % surface trace [pos dir] = surf_trace(pos_in,dir_in,cv,th,rno,rni); % draw surface profile y = linspace(-13,13,101); z = sag(y,cv); plot(z,y,'k','LineWidth',1.5); axis equal hold on; % extend ray three units after refraction ray = pos_in; ray(2,:)=pos; ray(3,:)=pos+3*dir; % draw rays plot(ray(:,3),ray(:,2),'b'); % surface normal vector ns = -cv*pos; ns(3)=1+ns(3); p = pos; p(2,:)=pos-5*ns; plot(p(:,3),p(:,2),'g'); hold off;
Rayfan
cv = 1/34.53; % surface curvature rno=1.670028; % index of refraction after surface (out) rni=1; % index of refraction before surface (in) th=0; % distance from last surface pos_in = [0 3 -10]; % starting position of ray % draw surface profile y = linspace(-13,13,101); z = sag(y,cv); plot(z,y,'k','LineWidth',1.5); axis equal hold on; ray = pos_in; for ky=linspace(-0.6,0.6,11) dir_in = [0,ky,1]; dir_in = dir_in/norm(dir_in); [pos dir] = surf_trace(pos_in,dir_in,cv,th,rno,rni); ray(2,:)=pos; ray(3,:)=pos+3*dir; plot(ray(:,3),ray(:,2),'b'); end hold off;
Example
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 ];
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);
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 surface
fprintf('image plane: (y = %.6f, x = %.6f)\n',pos(2),pos(1));
image plane: (y = -0.004580, x = 0.000000)
OPD Example
Add a non-refraction reference surface as the last surface in the optical system.
% define 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));
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)
The results are
disp(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);
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
disp(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.5875562e-3; fprintf('total opd %g wave\n',opd/wvl);
opd =
0 -2.3419 9.8366 -5.7527 -1.7356
total opd 0.00629975
total opd 10.722 wave