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;