theta = linspace(0,pi/2);

n = 1.5;



ddelta = abs(sin(theta).*cos(theta)./sqrt(n^2-sin(theta).^2));

plot(180/pi*theta,ddelta);