i=1;
for theta_i=0 : 0.01 : 90

theta_t=asind(sin(deg2rad(theta_i))/1.5);
r_s=sin(deg2rad(theta_t-theta_i))/sin(deg2rad(theta_t+theta_i));

theta(i)=theta_i;
r(i)=r_s;
i=i+1;
end
plot(theta,r);