1 files changed, 79 insertions, 0 deletions
diff --git a/coupling_angles.m b/coupling_angles.m
new file mode 100644
index 0000000..3af9c87
--- /dev/null
+++ b/coupling_angles.m
@@ -0,0 +1,79 @@
+%% 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');