function Y = trig_ray(Qin,Uin,cv,th,rn)
%
%  TRIG_RAY  trigonometric raytrace
%
m = length(th);

Q = Qin;
U = Uin;

for (n=2:m)
   I = asin(sin(U)+Q*cv(n));
   Iprime = asin(rn(n-1)*sin(I)/rn(n));
   Uprime = U + Iprime - I;
   G = Q./(cos(U)+cos(I));
   Qprime = G.*(cos(Uprime)+cos(Iprime));
   beta = Iprime - Uprime;
   Y = G.*(1+cos(beta));
   % uncomment following for intermediate results
   % [n Y U I Iprime Q Qprime G]
   Q = Qprime + th(n)*sin(Uprime);
   U = Uprime;
end