%% Calculates incident angle for proper coupling into the disc via prism

% angle of the prism faces in degrees
prism_angle=45; 

%% prism index of refraction
% Rutile (TiO2) see http://refractiveindex.info/?group=CRYSTALS&material=TiO2
n_rutile_o = 2.4885;
n_rutile_e = 2.75324;
n_p=n_rutile_o

%% disk material index of refraction
% Magnesium Fluoride (MgF2)  see http://refractiveindex.info/?group=CRYSTALS&material=MgF2
n_MgF2_o = 1.3751;
n_MgF2_e = 1.38679;

n_d=n_MgF2_o

%% critical angle for beam from prism to disk
% recall n_d*sin(theta_d)=n_p*sin(theta_p) where angles are counted from normal to the face
% and we want theta_d to be 90 degrees for the total internal reflection
theta_p=asin(n_d/n_p);
% convert to degrees
theta_p_in_degrees=theta_p*180/pi

%% now lets see what angle it does with other face of the prism
theta_p_in_degrees_2=-(theta_p_in_degrees-prism_angle)

%% now we calculate refracted angle out of the prism into the air
theta_air=asin(n_p*sin(theta_p_in_degrees_2/180*pi));
% convert to degrees
theta_air_in_degrees=theta_air*180/pi


%% Lets make a picture
% 1st face of the prism lays at y=0 and spans from x=-1 to x=1
x_face_1=linspace(-1,1);
y_face_1=0*x_face_1;
% 2nd face to the right of origin  and at angle prism_angle with respect to negative x direction
x_face_2=linspace(1,0);
y_face_2=(1-x_face_2)*tan(prism_angle/180*pi);
% 3rd face to the left of origin  and at angle prism_angle with respect to negative x direction
x_face_3=linspace(-1,0);
y_face_3=(x_face_3+1)*tan(prism_angle/180*pi);

%% draw prism
figure(1); hold off;
plot(x_face_1, y_face_1, 'k-', x_face_2, y_face_2, 'k-', x_face_3, y_face_3, 'k-');
hold on

%% disk center will be located in point (0,-R);
R=.5
dc_x=0; dc_y=-R;
x_disk=dc_x+R*cos(linspace(0,2*pi));
y_disk=dc_y+R*sin(linspace(0,2*pi));
plot(x_disk,y_disk,'k-')


%% beam trace inside the prism
% crossing of the beam with 1st (right) face of the prism
% we are solving cot(theta_p)*x=cot(prism_angle)*(1-x)
x_cross_1=cot(prism_angle/180*pi)/(cot(prism_angle/180*pi)+cot(theta_p));
x_beam_prism_1=linspace(0,x_cross_1);
y_beam_prism_1=linspace(0,x_cross_1)*cot(theta_p);
plot(x_beam_prism_1, y_beam_prism_1, 'r-');

%% beam out of prism
x_out_strt=x_cross_1;
y_out_strt=x_cross_1*cot(theta_p);

% angle in the air relative to horizon
theta_air_rlh=(theta_air_in_degrees-prism_angle)/180*pi;
if (abs(theta_air_rlh)<pi/2) x_out_stop=1; else x_out_stop=0; end
x_out_1=linspace(x_out_strt, x_out_stop);
y_out_1=y_out_strt+tan(theta_air_rlh)*(x_out_1-x_out_strt);
plot(x_out_1, y_out_1, 'r-');

% force same aspect ratio for axis
axis('square');