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);