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

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

prism_angle = prism_angle_in_degrees*pi/180; 

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

n_prism=n_rutile_o ; % for horizontal or p polarization

%% 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;

% Measured 
n_MgF2_p = 1.465; % p-polarization


coupling_description='Rutile prism, MgF_{2} disk, p-polarization';
n_disk=n_MgF2_p;

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

%% now lets see what angle it does with other face of the prism
theta_prism_2=(prism_angle - theta_prism);
theta_prism_in_degrees_2=theta_prism_2*180/pi

%% now we calculate refracted angle out of the prism into the air with respect to the normal
% positive means above the normal
asin_arg=n_prism*sin(theta_prism_2);
if (abs(asin_arg)>1) error('quiting: at the right prism face  we experienced total internal reflection'); end
theta_air=asin(asin_arg);
% convert to degrees
theta_air_in_degrees=theta_air*180/pi

%% angle in the air relative to horizon
theta_air_rlh=(theta_air+prism_angle);
theta_air_rlh_in_degrees=theta_air_rlh*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);
% 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);

%% 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=.25;
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_prism)*x=tan(prism_angle)*(1-x)
x_cross_1=tan(prism_angle)/(tan(prism_angle)+cot(theta_prism));
y_cross_1=x_cross_1*cot(theta_prism);
x_beam_prism_1=linspace(0,x_cross_1);
y_beam_prism_1=x_beam_prism_1*cot(theta_prism);
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_prism);

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-');

%% draw normal at the point of 1st face intersection
theta_norm_rlh=pi/2-prism_angle;
% coordinates of normal outside the prism
x_norm_out=linspace(x_cross_1,1);
y_norm_out=y_cross_1+(x_norm_out-x_cross_1)*tan(theta_norm_rlh);

% coordinates of normal inside the prism
x_norm_in=linspace(x_cross_1,0);
y_norm_in=y_cross_1+(x_norm_in-x_cross_1)*tan(theta_norm_rlh);
plot(x_norm_out, y_norm_out, 'b-', x_norm_in, y_norm_in, 'b-');

%% arc to show theta_air
R_arc=.1; 
phi_arc=linspace(theta_air_rlh,theta_norm_rlh);
x_arc_air=x_cross_1+R_arc*cos(phi_arc);
y_arc_air=y_cross_1+R_arc*sin(phi_arc);
plot(x_arc_air, y_arc_air, 'b-');
%% annotation theta_air_in_degrees
str=sprintf('theta_{a}=%.1f',theta_air_in_degrees);
str=strcat('\',str,'^{o}');
text(x_arc_air(end)+.05,y_arc_air(end), str);

%% normal at the disk contact
y_norm_in=linspace(0.,0.5);
x_norm_in=y_norm_in*0;
plot(x_norm_in, y_norm_in, 'b-');

%% arc to show theta_prism
R_arc=.1; 
phi_arc=linspace(pi/2-theta_prism,pi/2);
x_arc_prism=R_arc*cos(phi_arc);
y_arc_prism=R_arc*sin(phi_arc);
plot(x_arc_prism, y_arc_prism, 'b-');
%% annotation theta_air_in_degrees
str=sprintf('theta_{p}=%.1f',theta_prism_in_degrees);
str=strcat('\',str,'^{o}');
text(x_arc_prism(1)+0.05, y_arc_prism(1), str);

%% arc to show theta_prism_2
R_arc=.1; 
phi_arc=linspace(theta_norm_rlh+theta_prism_2+pi,theta_norm_rlh+pi);
x_arc_prism_2=x_cross_1+R_arc*cos(phi_arc);
y_arc_prism_2=y_cross_1+R_arc*sin(phi_arc);
plot(x_arc_prism_2, y_arc_prism_2, 'b-');
%% annotation theta_air_in_degrees
str=sprintf('theta_{p2}=%.1f',theta_prism_in_degrees_2);
str=strcat('\',str,'^{o}');
text(x_arc_prism_2(1)+.05,y_arc_prism_2(1), str);

%% arc to show prism angle
R_arc=.1; 
phi_arc=linspace(pi, pi-prism_angle);
x_arc=1+R_arc*cos(phi_arc);
y_arc=0+R_arc*sin(phi_arc);
plot(x_arc, y_arc, 'b-');
%% annotate prism angle
str=sprintf('theta_{prism}=%.1f',prism_angle_in_degrees);
str=strcat('\',str,'^{o}');
text(.65, 0.06, str);


%% General annotation
str=sprintf('n_{disk}=%.3f',n_disk);
text(-0.9,1.1, str);

str=sprintf('n_{prism}=%.3f',n_prism);
text(-0.9,1.2, str);

text(-0.9, 1.3, coupling_description);

% force same aspect ratio for axis
axis([-1,1,-0.5,1.5],'equal');

%% output of the plot to the file
print('prism_disk_coupling.eps','-depsc2');
print('prism_disk_coupling.png');