parax

function pray = parax(pstart, cv, th, rn)
%
%PARAX
%
% paraxial ray trace

pray(1,:) = pstart(1,:);
m = length(th);

if (nargin==3)
   for (n=2:m)
      pray(n,1) = pray(n-1,1)+ pray(n-1,2)*th(n-1);
      power = cv(n);
      pray(n,2) = pray(n-1,2) - power*pray(n,1);
   end
else
   nu = rn(1)*pray(1,2);
   for (n=2:m)
      pray(n,1) = pray(n-1,1)+ pray(n-1,2)*th(n-1);
      power = (rn(n)-rn(n-1))*cv(n);
      nu = nu - power * pray(n,1);
      pray(n,2) = nu / rn(n);
   end
end
pray(m+1,1) = pray(m,1)+th(m)*pray(m,2);
pray(m+1,2) = pray(m,2);