summaryrefslogtreecommitdiff
path: root/beam_tracing/octave
diff options
context:
space:
mode:
Diffstat (limited to 'beam_tracing/octave')
-rw-r--r--beam_tracing/octave/README10
-rw-r--r--beam_tracing/octave/TODO3
-rw-r--r--beam_tracing/octave/beam2face_distance.m47
-rw-r--r--beam_tracing/octave/beam_trace.m47
-rw-r--r--beam_tracing/octave/color_gradient.m31
-rw-r--r--beam_tracing/octave/demo_prism_disk_coupling.m29
-rw-r--r--beam_tracing/octave/example_beam_trace.m75
-rw-r--r--beam_tracing/octave/example_beam_trace_diamond.m76
-rw-r--r--beam_tracing/octave/face_beam_interaction.m151
-rw-r--r--beam_tracing/octave/fresnel_reflection.m28
-rw-r--r--beam_tracing/octave/make_beam_trace.m37
-rw-r--r--beam_tracing/octave/plot_beams_and_faces_figure.m56
-rw-r--r--beam_tracing/octave/prism_disk_coupling.m152
13 files changed, 742 insertions, 0 deletions
diff --git a/beam_tracing/octave/README b/beam_tracing/octave/README
new file mode 100644
index 0000000..4fe8125
--- /dev/null
+++ b/beam_tracing/octave/README
@@ -0,0 +1,10 @@
+= WGMR related code =
+
+prism_disk_coupling.m - calculates parameters of the beam to be coupled
+ into the disk via prism
+demo_prism_disk_coupling.m - demo case of the above routine
+beam_trace.m - calculates beams propagation and their branching
+ on dielectric interfaces (faces) for arbitrary set of beams and faces.
+example_beam_trace.m - example use of above routine.
+
+rest of the code: helper functions
diff --git a/beam_tracing/octave/TODO b/beam_tracing/octave/TODO
new file mode 100644
index 0000000..8265edc
--- /dev/null
+++ b/beam_tracing/octave/TODO
@@ -0,0 +1,3 @@
+* Add an plotting functions which plot only beams and faces
+ after calculations of all beams done.
+
diff --git a/beam_tracing/octave/beam2face_distance.m b/beam_tracing/octave/beam2face_distance.m
new file mode 100644
index 0000000..d7d5a61
--- /dev/null
+++ b/beam_tracing/octave/beam2face_distance.m
@@ -0,0 +1,47 @@
+function [hit_distance, hit_position, is_face_hit] = beam2face_distance(beam,face)
+ %% return distance to face if beam hits it or infinity otherwise
+ % for beams and faces structure description see face_beam_interaction.m
+
+ % default returning values for the case when beam misses the face
+ is_face_hit = false;
+ hit_position = [NA, NA];
+ hit_distance = Inf;
+
+ k = beam.k;
+ % face direction or kf
+ kf=face.vertex2 - face.vertex1; % not a unit vector
+
+ %% simple check for intersection of two vectors
+ % if beam are parallel no intersection is possible
+ % we do this via simulated cross product calculation
+ angle_tolerance = 1e-15;
+ if ( abs( (k(1)*kf(2)-k(2)*kf(1))/norm(k)/norm(kf) ) <= angle_tolerance )
+ % beams never intercept
+ return;
+ end
+
+ %% let's find the intersection of lines passing through face and light beam
+ % we introduce parameters tb and tf
+ % which define distance (in arb. units) along k vectors
+ % we are solving
+ % xb_o+ k_x*tb = v1_x+ kf_x*tf
+ % yb_o+ k_y*tb = v1_y+ kf_y*tf
+ % A*t = B
+ A = [ k(1) , -kf(1); k(2), -kf(2)];
+ B= [ face.vertex1(1)-beam.origin(1); face.vertex1(2)-beam.origin(2) ];
+ t = A\B; tb=t(1); tf=t(2);
+
+ % beam intersects face only if 0<=tf<=1 and tb>=0
+ if ( (0 <= tf) && (tf<=1) && tb >=0 )
+ % intersection, beam hits the face
+ is_face_hit = true;
+ else
+ % beam misses the face
+ return;
+ end
+
+ % calculating hit position
+ hit_position=face.vertex1+kf*tf;
+ hit_distance = tb; % distance along the light beam
+end
+
diff --git a/beam_tracing/octave/beam_trace.m b/beam_tracing/octave/beam_trace.m
new file mode 100644
index 0000000..90359d4
--- /dev/null
+++ b/beam_tracing/octave/beam_trace.m
@@ -0,0 +1,47 @@
+function [processed_beams] = beam_trace(beams, faces, borders, border_limits)
+ %% trace beam and all it reflections and refractions on faces
+ % for beams and faces structure description see face_beam_interaction.m
+ % border similar to faces cell array which enclose the beam
+ % i.e. it does not leave the polygon consisting of border faces
+
+ intensity_threshold = 1e-3; % intensity of the weakest beam which we still trace
+ processed_beams = {};
+
+ Nbeams=size(beams)(2);
+ if ( Nbeams == 0 )
+ % no more beam to trace
+ return;
+ end
+
+ % trace each beam in beams cell array
+ for i=1:Nbeams
+ beam=beams{i};
+ [is_face_hit, hit_position, hit_distance, new_beams] = face_beam_interaction(beam, faces);
+ if (is_face_hit)
+ beams{i}.hit_position=hit_position;
+ %% remove beams which are to weak to trace
+ N_new_beams=size(new_beams)(2);
+ intense_enough_beams={};
+ intense_beams_counter=0;
+ for k=1:N_new_beams
+ if ( new_beams{k}.intensity >= intensity_threshold )
+ intense_beams_counter=intense_beams_counter + 1;
+ intense_enough_beams{intense_beams_counter}=new_beams{k};
+ end
+ end
+ [new_processed_beams] = beam_trace(intense_enough_beams, faces, borders, border_limits);
+ processed_beams=cat(2,processed_beams, new_processed_beams);
+ else
+ % beam does not hit face but it should stop and borders
+ beam.face=NA;
+ [is_face_hit, hit_position, hit_distance, new_beams] = face_beam_interaction(beam, borders);
+ if (!is_face_hit)
+ error('borders are badly defined, the beam misses them');
+ end
+ beams{i}.hit_position=hit_position;
+ end
+ end
+ processed_beams=cat(2,beams, processed_beams);
+end
+
+
diff --git a/beam_tracing/octave/color_gradient.m b/beam_tracing/octave/color_gradient.m
new file mode 100644
index 0000000..310ff99
--- /dev/null
+++ b/beam_tracing/octave/color_gradient.m
@@ -0,0 +1,31 @@
+function rgb= color_gradient( pos, base_color_rgb, Ncolors )
+%% returns rgb=[R,G,B] components from the brightness gradient of the base color
+% base_color_rgb - [R,G,B] of the base color
+% Ncolors - number of positions in gradient
+
+ %% color map as gradient of particular color
+ % starting values
+ hsv_s=rgb2hsv(base_color_rgb);
+ hue_f = hsv_s(1);
+ sat_f = hsv_s(2);
+ val_f = hsv_s(3);
+ % stop values
+ hue_s = hue_f;
+ sat_s = 0;
+ val_s = 1;
+ hue=linspace(hue_s,hue_f, Ncolors);
+ sat=linspace(sat_s,sat_f, Ncolors);
+ val=linspace(val_s,val_f, Ncolors);
+ cmap = hsv2rgb([hue', sat', val']);
+
+ %% check that pos is with in the reach
+ pos = floor(pos);
+ if( pos < 1 )
+ pos=1;
+ end
+ if( pos > Ncolors )
+ pos=Ncolors;
+ end
+
+ rgb = cmap(pos,:);
+end
diff --git a/beam_tracing/octave/demo_prism_disk_coupling.m b/beam_tracing/octave/demo_prism_disk_coupling.m
new file mode 100644
index 0000000..e809968
--- /dev/null
+++ b/beam_tracing/octave/demo_prism_disk_coupling.m
@@ -0,0 +1,29 @@
+%% 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
+
+
+n_disk=n_MgF2_p;
+coupling_description='Rutile prism, MgF_{2} disk, p-polarization';
+
+% calculate coupling angle and draw coupling solution with annotations
+prism_disk_coupling(prism_angle_in_degrees, n_disk, n_prism, coupling_description);
+
diff --git a/beam_tracing/octave/example_beam_trace.m b/beam_tracing/octave/example_beam_trace.m
new file mode 100644
index 0000000..6827e54
--- /dev/null
+++ b/beam_tracing/octave/example_beam_trace.m
@@ -0,0 +1,75 @@
+%beam1.k=[-1,0];
+%beam1.origin=[2,0.5];
+beam_initial_travel_angle=(28+45-180)/180*pi;
+beam1.k=[cos(beam_initial_travel_angle), sin(beam_initial_travel_angle)];
+beam1.origin=[1.0,2];
+%beam1.k=[-1,-1];
+%beam1.origin=[2,2];
+beam1.face=NA;
+beam1.intensity = 1;
+beam1.polarization = 1; % 1 for s and 2 for p
+beam1.status='incoming'; % could be reflected, refracted, incoming
+
+beam2=beam1;
+beam1.polarization = 2; % 1 for s and 2 for p
+
+beams={beam1,beam2};
+
+n_rutile=[2.7532, 2.4885]; % n_s and n_p
+n_disk= [1.3751, 1.3868]; % MgF2 book value
+n_disk= [1.465 , 1.465 ]; % MgF2 measured
+n_disk= [1.39 , 1.39 ]; % CaF2 measured for p-polarization
+n_air= [1., 1.];
+
+face1.vertex1=[-1,0];
+face1.vertex2=[1,0];
+%face1.n_right = n_air;
+face1.n_right = n_disk;
+face1.n_left= n_rutile;
+
+face2.vertex1=[1,0];
+face2.vertex2=[0,1];
+face2.n_right = [1, 1];
+face2.n_left= n_rutile;
+
+face3.vertex1=[0,1];
+face3.vertex2=[-1,0];
+face3.n_right = [1, 1];
+face3.n_left= n_rutile;
+
+faces={face1,face2,face3};
+
+
+border1.vertex1=[-2,-2];
+border1.vertex2=[2,-2];
+border1.n_right = [1, 1];
+border1.n_left= [1, 1];
+
+border2.vertex1=[2,-2];
+border2.vertex2=[2,2];
+border2.n_right = [1, 1];
+border2.n_left= [1, 1];
+
+border3.vertex1=[2,2];
+border3.vertex2=[-2,2];
+border3.n_right = [1, 1];
+border3.n_left= [1, 1];
+
+border4.vertex1=[-2,2];
+border4.vertex2=[-2,-2];
+border4.n_right = [1, 1];
+border4.n_left= [1, 1];
+
+borders={border1,border2, border3, border4};
+
+
+border_limits=[-2,-2, 2,2];
+figure(1);
+clf;
+[processed_beams] = beam_trace(beams, faces, borders, border_limits );
+text(-1,-1, 'red: s-polarization; blue: p-polarization');
+
+
+print('beams.eps','-depsc2');
+print('beams.png');
+
diff --git a/beam_tracing/octave/example_beam_trace_diamond.m b/beam_tracing/octave/example_beam_trace_diamond.m
new file mode 100644
index 0000000..e3cc30d
--- /dev/null
+++ b/beam_tracing/octave/example_beam_trace_diamond.m
@@ -0,0 +1,76 @@
+%beam1.k=[-1,0];
+%beam1.origin=[2,0.5];
+beam_initial_travel_angle=(70+30-180)/180*pi;
+beam1.k=[cos(beam_initial_travel_angle), sin(beam_initial_travel_angle)];
+beam1.origin=[0.5,2];
+%beam1.k=[-1,-1];
+%beam1.origin=[2,2];
+beam1.face=NA;
+beam1.intensity = 1;
+beam1.polarization = 1; % 1 for s and 2 for p
+beam1.status='incoming'; % could be reflected, refracted, incoming
+
+beam2=beam1;
+beam1.polarization = 2; % 1 for s and 2 for p
+
+beams={beam1,beam2};
+
+n_rutile=[2.7532, 2.4885]; % n_s and n_p
+n_diamond=[2.4 , 2.4 ]; % n_s and n_p
+n_disk= [1.3751, 1.3868]; % MgF2 book value
+n_disk= [1.465 , 1.465 ]; % MgF2 measured
+n_disk= [1.39 , 1.39 ]; % CaF2 measured for p-polarization
+n_air= [1., 1.];
+
+face1.vertex1=[-1,0];
+face1.vertex2=[1,0];
+%face1.n_right = n_air;
+face1.n_right = n_disk;
+face1.n_left= n_diamond;
+
+face2.vertex1=[1,0];
+face2.vertex2=[0,sqrt(3)];
+face2.n_right = [1, 1];
+face2.n_left= n_diamond;
+
+face3.vertex1=[0,sqrt(3)];
+face3.vertex2=[-1,0];
+face3.n_right = [1, 1];
+face3.n_left= n_diamond;
+
+faces={face1,face2,face3};
+
+
+border1.vertex1=[-2,-2];
+border1.vertex2=[2,-2];
+border1.n_right = [1, 1];
+border1.n_left= [1, 1];
+
+border2.vertex1=[2,-2];
+border2.vertex2=[2,2];
+border2.n_right = [1, 1];
+border2.n_left= [1, 1];
+
+border3.vertex1=[2,2];
+border3.vertex2=[-2,2];
+border3.n_right = [1, 1];
+border3.n_left= [1, 1];
+
+border4.vertex1=[-2,2];
+border4.vertex2=[-2,-2];
+border4.n_right = [1, 1];
+border4.n_left= [1, 1];
+
+borders={border1,border2, border3, border4};
+
+
+border_limits=[-2,-2, 2,2];
+figure(1);
+clf;
+[processed_beams] = beam_trace(beams, faces, borders, border_limits );
+text(-1,-1, 'red: s-polarization; blue: p-polarization');
+
+
+print('beams.eps','-depsc2');
+print('beams.png');
+
diff --git a/beam_tracing/octave/face_beam_interaction.m b/beam_tracing/octave/face_beam_interaction.m
new file mode 100644
index 0000000..3bc3e4a
--- /dev/null
+++ b/beam_tracing/octave/face_beam_interaction.m
@@ -0,0 +1,151 @@
+function [is_face_hit, hit_position, hit_distance, new_beams] = face_beam_interaction(beam, faces)
+ %% calculates refracted and reflected beam after interaction with a face
+ % beam - structure defining the light beam
+ % beam.origin - array [x,y] origin/soutce of the light beam
+ % beam.k - k vector i.e. direction [kx,ky]
+ % beam.intensity - intensity of the beam
+ % beam.face - if beam starts from face then its index is here
+ % beam.polarization - can be either 1 for 's' or 2 for 'p'
+ % beam.status - could be 'reflected', 'refracted', 'incoming'
+ % faces cell array of face structures
+ % face - structure definiong the beam
+ % face.vertex1 - [x,y] of the 1st point/vertex of the face
+ % face.vertex2 - [x,y] of the 2nd point/vertex of the face
+ % face.n_left - indexes of refraction on the left hand side
+ % with respect to 1st -> 2nd vertex direction
+ % [ns,np] - for s and p polarizations
+ % face.n_right - indexes of refraction on the right hand side
+ % [ns,np] - for s and p polarizations
+
+
+
+ k=beam.k;
+ polarization=beam.polarization;
+
+ %% we go over all faces to find the closest which beam hits
+ Nfaces=size(faces)(2);
+ hit_distance=Inf;
+ is_face_hit = false;
+ hit_position = [NA, NA];
+ closest_face_index=NA;
+ for i=1:Nfaces
+ if ( beam.face == i) continue; end
+ face=faces{i};
+ [hit_distance_tmp, hit_position_tmp, is_face_hit_tmp] = beam2face_distance(beam,face);
+ if ( hit_distance_tmp < hit_distance )
+ % this is the closer face
+ is_face_hit=is_face_hit_tmp;
+ hit_position=hit_position_tmp;
+ hit_distance=hit_distance_tmp;
+ closest_face_index=i;
+ end
+ end
+
+ if (!is_face_hit)
+ new_beams={};
+ return;
+ end
+
+ %% closest face
+ face=faces{closest_face_index};
+ kf=face.vertex2 - face.vertex1; % not a unit vector
+
+ hold on;
+ % draw face
+ line([face.vertex1(1), face.vertex2(1)], [face.vertex1(2), face.vertex2(2)] , 'linestyle', '-', 'color', [0,0,0] );
+ % draw beam
+ Ncolors=256;
+ linewidth=2;
+ if ( polarization == 1)
+ % s-polarization
+ line_base_color=[1,0,0]; % RGB - red
+ else
+ % p-polarization
+ line_base_color=[0,0,1]; % RGB - blue
+ end
+ %plot(x,y,line_str);
+ line_color = color_gradient( beam.intensity*Ncolors, line_base_color, Ncolors);
+ line([beam.origin(1), hit_position(1)], [beam.origin(2), hit_position(2)], 'color', line_color, 'linewidth', 3, 'linestyle', '-');
+
+ % find is beam arriving from left or right. I will use vectors cross product property.
+ % if z component of the product 'k x kf' is positive then beam arrives from the left
+ if ( ( k(1)*kf(2)-k(2)*kf(1) ) > 0 )
+ % beam coming from the left
+ n1=face.n_left(polarization);
+ n2=face.n_right(polarization);
+ % this means that notmal and beam k vector point to the same half plane
+ % relative to the face
+ are_nf_and_k_in_the_same_plane=true;
+ else
+ % beam coming from the right
+ n1=face.n_right(polarization);
+ n2=face.n_left(polarization);
+ % this means that notmal and beam k vector point to the different half plane
+ % relative to the face
+ are_nf_and_k_in_the_same_plane=false;
+ end
+
+ % normal vector to the face, looks to the left of it
+ nf=[ kf(2), -kf(1) ] / norm(kf);
+ % incidence angle calculation
+ cos_theta_i = dot(k, nf) / (norm(k)*norm(nf));
+ sin_theta_i = - ( k(1)*nf(2)-k(2)*nf(1) ) / (norm(k)*norm(nf));
+ % positive angle to the right from normal before incidence to the face
+ theta_i = atan2(sin_theta_i, cos_theta_i);
+
+ % we need to make sure that angle of incidence belong to [-pi/2,pi/2] interval
+ if( theta_i > pi/2)
+ theta_i=pi-theta_i;
+ end
+ if( theta_i < -pi/2)
+ theta_i=-pi-theta_i;
+ end
+
+ % angle of the normal with respect to horizon
+ theta_normal = atan2(nf(2), nf(1));
+
+ %% reflected beam direction
+ if (are_nf_and_k_in_the_same_plane)
+ theta_reflected = theta_normal - theta_i + pi;
+ else
+ theta_reflected = theta_normal + theta_i;
+ end
+
+ %% coefficients of reflection and transmission for given polarization
+ [R,theta_refracted_rel2normal]=fresnel_reflection(n1, n2, theta_i);
+ reflectivity = R(polarization);
+ transmission = 1 - reflectivity;
+
+ beam_reflected.origin = hit_position;
+ beam_reflected.k = [cos(theta_reflected), sin(theta_reflected)];
+ beam_reflected.face=closest_face_index;
+ beam_reflected.polarization = polarization;
+ beam_reflected.intensity = beam.intensity * reflectivity;
+ beam_reflected.status = 'reflected';
+ new_beams{1} = beam_reflected;
+
+
+ %% refracted beam direction
+ % refracted angle with respect to normal
+ if ( transmission == 0 )
+ % total internal reflection
+ else
+ % beam refracts
+ if (are_nf_and_k_in_the_same_plane)
+ theta_refracted = theta_normal + theta_refracted_rel2normal;
+ else
+ theta_refracted = theta_normal - theta_refracted_rel2normal + pi;
+ end
+
+ beam_refracted.origin = hit_position;
+ beam_refracted.k = [cos(theta_refracted), sin(theta_refracted)];
+ beam_refracted.face=closest_face_index;
+ beam_refracted.polarization = polarization;
+ beam_refracted.intensity = beam.intensity * transmission;
+ beam_refracted.status = 'refracted';
+ new_beams{2} = beam_refracted;
+ end
+end
+
+
+
diff --git a/beam_tracing/octave/fresnel_reflection.m b/beam_tracing/octave/fresnel_reflection.m
new file mode 100644
index 0000000..347d60b
--- /dev/null
+++ b/beam_tracing/octave/fresnel_reflection.m
@@ -0,0 +1,28 @@
+function [R, theta_t] = fresnel_reflection(n1, n2, theta_i)
+%% calculates intensity reflection coefficient for s and p polarizations
+%% for light travelling from material with index of refraction n1 to material with n2
+%% theta_i - incident angle in medium 1 with respect to normal, could be a vector
+%% R - coefficients of reflection array [Rs, Rp]
+%% theta_t - transmitted/refracted angle in medium 2 with respect to normal
+
+ if ( size(theta_i)(1) != 1)
+ error('theta_i must be a vector or scalar');
+ end
+
+ %% see http://en.wikipedia.org/wiki/Fresnel_equations
+ %% refraction angle or angle of transmitted beam with respect to normal
+ sin_theta_t=n1/n2*sin(theta_i);
+ %% special cases: total internal reflection
+ indx = (abs(sin_theta_t) >= 1);
+ sin_theta_t( indx ) = sign( sin_theta_t(indx) );
+
+
+ theta_t = asin(sin_theta_t); % angle of refraction or transmitted beam
+
+ cos_theta_t = cos( theta_t );
+ cos_theta_i = cos( theta_i );
+
+ Rs = ( ( n1*cos_theta_i - n2*cos_theta_t ) ./ ( n1*cos_theta_i + n2*cos_theta_t ) ).^2;
+ Rp = ( ( n1*cos_theta_t - n2*cos_theta_i ) ./ ( n1*cos_theta_t + n2*cos_theta_i ) ).^2;
+ R = [Rs; Rp];
+end
diff --git a/beam_tracing/octave/make_beam_trace.m b/beam_tracing/octave/make_beam_trace.m
new file mode 100644
index 0000000..f62168f
--- /dev/null
+++ b/beam_tracing/octave/make_beam_trace.m
@@ -0,0 +1,37 @@
+function img=make_beam_trace(beam, stop_point, border_limits, img)
+ % img so far beam traced part
+ [Ny,Nx]=size(img);
+ %% border_limits has coordinates of left bottom and right top coners
+ xlb=border_limits(1);
+ ylb=border_limits(2);
+ xrt=border_limits(3);
+ yrt=border_limits(4);
+
+ %% beam start stop coordinates with respect to lower left border point
+ xb1=beam.origin(1)-xlb;
+ yb1=beam.origin(2)-ylb;
+ xb2=stop_point(1)-xlb;
+ yb2=stop_point(2)-ylb;
+
+ %% beam path coordinates
+ Nc=1000; % number of coordinate points
+ xb=linspace(xb1,xb2, Nc);
+ yb=linspace(yb1,yb2, Nc);
+
+ %% beam coordinate to image pixel coordinate
+ px=round(xb/(xrt-xlb)*(Nx-1) + 1);
+ py=round(yb/(yrt-ylb)*(Ny-1) + 1);
+
+ %% truncate to borders
+ py=py( (1<=px) & (px<=Nx) );
+ px=px( (1<=px) & (px<=Nx) );
+
+ px=px( (1<=py) & (py<=Ny) );
+ py=py( (1<=py) & (py<=Ny) );
+
+ for i=1:length(px)
+ img(py(i),px(i))=beam.intensity;
+ end
+end
+
+
diff --git a/beam_tracing/octave/plot_beams_and_faces_figure.m b/beam_tracing/octave/plot_beams_and_faces_figure.m
new file mode 100644
index 0000000..a63ac85
--- /dev/null
+++ b/beam_tracing/octave/plot_beams_and_faces_figure.m
@@ -0,0 +1,56 @@
+function plot_beams_and_faces_figure(border_limits, img, faces, fig_handle)
+ %% plot faces and beams images
+ % border_limits has coordinates of left bottom and right top coners
+ xlb=border_limits(1);
+ ylb=border_limits(2);
+ xrt=border_limits(3);
+ yrt=border_limits(4);
+
+ [Ny,Nx]=size(img);
+ xc=linspace(xlb,xrt, Nx);
+ yc=linspace(ylb,yrt, Ny);
+
+ figure(fig_handle);
+ %% dummy plot just to put axis in the proper directions
+ plot(xlb,ylb,'.', xrt,yrt, '.');
+ hold on;
+
+
+ %% color map as gradient of particular color
+ Ncolors=64; % number of colors
+ % starting values
+ hsv_s=rgb2hsv([1,0,0]);
+ hue_f = hsv_s(1);
+ sat_f = hsv_s(2);
+ val_f = hsv_s(3);
+ % stop values
+ hue_s = hue_f;
+ sat_s = 0;
+ val_s = 1;
+ hue=linspace(hue_s,hue_f, Ncolors);
+ sat=linspace(sat_s,sat_f, Ncolors);
+ val=linspace(val_s,val_f, Ncolors);
+ cmap = hsv2rgb([hue', sat', val']);
+
+ colormap( cmap );
+
+ %% plot ray images
+ imagesc(xc,yc, img); colorbar;
+ hold on;
+
+ %% plot all faces
+ Nf=size(faces)(2); % number of faces
+ for i=1:Nf
+ hold on;
+ t=linspace(0,1,100);
+ xf=faces{i}.vertex1(1) + (faces{i}.vertex2(1)-faces{i}.vertex1(1))*t;
+ yf=faces{i}.vertex1(2) + (faces{i}.vertex2(2)-faces{i}.vertex1(2))*t;
+ plot(xf,yf, 'k-');
+ end
+
+ axis('equal');
+ hold off;
+
+
+
+end
diff --git a/beam_tracing/octave/prism_disk_coupling.m b/beam_tracing/octave/prism_disk_coupling.m
new file mode 100644
index 0000000..d89806c
--- /dev/null
+++ b/beam_tracing/octave/prism_disk_coupling.m
@@ -0,0 +1,152 @@
+function prism_disk_coupling(prism_angle_in_degrees, n_disk, n_prism, coupling_description)
+%% Calculates incident angle for proper coupling into the disc via prism
+% prism_angle_in_degrees - angle of the prism faces in degrees
+% coupling_description - short annotation of the situation
+% for example:
+% coupling_description='Rutile prism, MgF_{2} disk, p-polarization';
+
+prism_angle = prism_angle_in_degrees*pi/180;
+
+%% 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+( pi/2 - 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(.52, 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');
+
+end
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-16 08:33:18 +0000