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