From a08ed173692ce16b79002733003049ec32a02485 Mon Sep 17 00:00:00 2001 From: Eugeniy Mikhailov Date: Thu, 28 Feb 2013 14:58:35 -0500 Subject: moved octave using code to separate folder --- beam_tracing/README | 10 -- beam_tracing/TODO | 3 - beam_tracing/beam2face_distance.m | 47 ------- beam_tracing/beam_trace.m | 47 ------- beam_tracing/color_gradient.m | 31 ----- beam_tracing/demo_prism_disk_coupling.m | 29 ----- beam_tracing/example_beam_trace.m | 75 ----------- beam_tracing/example_beam_trace_diamond.m | 76 ----------- beam_tracing/face_beam_interaction.m | 151 --------------------- beam_tracing/fresnel_reflection.m | 28 ---- beam_tracing/make_beam_trace.m | 37 ------ beam_tracing/octave/README | 10 ++ beam_tracing/octave/TODO | 3 + beam_tracing/octave/beam2face_distance.m | 47 +++++++ beam_tracing/octave/beam_trace.m | 47 +++++++ beam_tracing/octave/color_gradient.m | 31 +++++ beam_tracing/octave/demo_prism_disk_coupling.m | 29 +++++ beam_tracing/octave/example_beam_trace.m | 75 +++++++++++ beam_tracing/octave/example_beam_trace_diamond.m | 76 +++++++++++ beam_tracing/octave/face_beam_interaction.m | 151 +++++++++++++++++++++ beam_tracing/octave/fresnel_reflection.m | 28 ++++ beam_tracing/octave/make_beam_trace.m | 37 ++++++ beam_tracing/octave/plot_beams_and_faces_figure.m | 56 ++++++++ beam_tracing/octave/prism_disk_coupling.m | 152 ++++++++++++++++++++++ beam_tracing/plot_beams_and_faces_figure.m | 56 -------- beam_tracing/prism_disk_coupling.m | 152 ---------------------- 26 files changed, 742 insertions(+), 742 deletions(-) delete mode 100644 beam_tracing/README delete mode 100644 beam_tracing/TODO delete mode 100644 beam_tracing/beam2face_distance.m delete mode 100644 beam_tracing/beam_trace.m delete mode 100644 beam_tracing/color_gradient.m delete mode 100644 beam_tracing/demo_prism_disk_coupling.m delete mode 100644 beam_tracing/example_beam_trace.m delete mode 100644 beam_tracing/example_beam_trace_diamond.m delete mode 100644 beam_tracing/face_beam_interaction.m delete mode 100644 beam_tracing/fresnel_reflection.m delete mode 100644 beam_tracing/make_beam_trace.m create mode 100644 beam_tracing/octave/README create mode 100644 beam_tracing/octave/TODO create mode 100644 beam_tracing/octave/beam2face_distance.m create mode 100644 beam_tracing/octave/beam_trace.m create mode 100644 beam_tracing/octave/color_gradient.m create mode 100644 beam_tracing/octave/demo_prism_disk_coupling.m create mode 100644 beam_tracing/octave/example_beam_trace.m create mode 100644 beam_tracing/octave/example_beam_trace_diamond.m create mode 100644 beam_tracing/octave/face_beam_interaction.m create mode 100644 beam_tracing/octave/fresnel_reflection.m create mode 100644 beam_tracing/octave/make_beam_trace.m create mode 100644 beam_tracing/octave/plot_beams_and_faces_figure.m create mode 100644 beam_tracing/octave/prism_disk_coupling.m delete mode 100644 beam_tracing/plot_beams_and_faces_figure.m delete mode 100644 beam_tracing/prism_disk_coupling.m diff --git a/beam_tracing/README b/beam_tracing/README deleted file mode 100644 index 4fe8125..0000000 --- a/beam_tracing/README +++ /dev/null @@ -1,10 +0,0 @@ -= 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/TODO b/beam_tracing/TODO deleted file mode 100644 index 8265edc..0000000 --- a/beam_tracing/TODO +++ /dev/null @@ -1,3 +0,0 @@ -* Add an plotting functions which plot only beams and faces - after calculations of all beams done. - diff --git a/beam_tracing/beam2face_distance.m b/beam_tracing/beam2face_distance.m deleted file mode 100644 index d7d5a61..0000000 --- a/beam_tracing/beam2face_distance.m +++ /dev/null @@ -1,47 +0,0 @@ -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/beam_trace.m b/beam_tracing/beam_trace.m deleted file mode 100644 index 90359d4..0000000 --- a/beam_tracing/beam_trace.m +++ /dev/null @@ -1,47 +0,0 @@ -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/color_gradient.m b/beam_tracing/color_gradient.m deleted file mode 100644 index 310ff99..0000000 --- a/beam_tracing/color_gradient.m +++ /dev/null @@ -1,31 +0,0 @@ -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/demo_prism_disk_coupling.m b/beam_tracing/demo_prism_disk_coupling.m deleted file mode 100644 index e809968..0000000 --- a/beam_tracing/demo_prism_disk_coupling.m +++ /dev/null @@ -1,29 +0,0 @@ -%% 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/example_beam_trace.m b/beam_tracing/example_beam_trace.m deleted file mode 100644 index 6827e54..0000000 --- a/beam_tracing/example_beam_trace.m +++ /dev/null @@ -1,75 +0,0 @@ -%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/example_beam_trace_diamond.m b/beam_tracing/example_beam_trace_diamond.m deleted file mode 100644 index e3cc30d..0000000 --- a/beam_tracing/example_beam_trace_diamond.m +++ /dev/null @@ -1,76 +0,0 @@ -%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/face_beam_interaction.m b/beam_tracing/face_beam_interaction.m deleted file mode 100644 index 3bc3e4a..0000000 --- a/beam_tracing/face_beam_interaction.m +++ /dev/null @@ -1,151 +0,0 @@ -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/fresnel_reflection.m b/beam_tracing/fresnel_reflection.m deleted file mode 100644 index 347d60b..0000000 --- a/beam_tracing/fresnel_reflection.m +++ /dev/null @@ -1,28 +0,0 @@ -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/make_beam_trace.m b/beam_tracing/make_beam_trace.m deleted file mode 100644 index f62168f..0000000 --- a/beam_tracing/make_beam_trace.m +++ /dev/null @@ -1,37 +0,0 @@ -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/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)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)