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
First we define the optical system.
>> cv = [ 0 0.034654 -0.005408 0]; >> th = [ 0 1 49.42301 0]; >> rn = [ 1 1.5 1 1 ];
Then we can generate a single trigonmetric ray:
>> format short g >> ye = trig_ray(6.25, 0, cv, th, rn); ye = -0.10864
We can use trig_ray to trace a fan of rays:
>> y = linspace(0,1,6)'; >> ye = trig_ray(y*6.25,0,cv,th,rn); >> [y ye] ans = 0 0 0.2 -0.00083418 0.4 -0.0067075 0.6 -0.02283 0.8 -0.054764 1 -0.10864
We can fit the wavefront polynomial as follows:
>> lambda = 0.00055; >> ua = -0.125; >> epsilon = abs(lambda/ua); >> d = - ye/epsilon; >> c = [ 2*y 4*y.^3 6*y.^5 ]; >> a = c\d a = 0.0022615 5.9053 0.17763
The rms deviation can be calculated as:
>> diff = d - c*a; >> rms = sqrt(sum(diff.*diff))/length(diff) rms = 0.00018571
Maintained by John Loomis, last updated 24 May 1999