From 2e6895ec0b70ba97e06c00f0447f86ad162181fd Mon Sep 17 00:00:00 2001 From: Eugeniy Mikhailov Date: Thu, 7 Jul 2011 01:21:13 -0400 Subject: arbitrary beam tracing --- beam2face_distance.m | 36 ++++++++++++++++ beam_trace.m | 31 ++++++++++++++ coupling_angles.m | 3 +- example.m | 30 ++++++++++++++ example_beam_trace.m | 45 ++++++++++++++++++++ face_beam_interaction.m | 107 ++++++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 251 insertions(+), 1 deletion(-) create mode 100644 beam2face_distance.m create mode 100644 beam_trace.m create mode 100644 example.m create mode 100644 example_beam_trace.m create mode 100644 face_beam_interaction.m diff --git a/beam2face_distance.m b/beam2face_distance.m new file mode 100644 index 0000000..ba3d16b --- /dev/null +++ b/beam2face_distance.m @@ -0,0 +1,36 @@ +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 + + k = beam.k; + % face direction or kf + kf=face.vertex2 - face.vertex1; % not a unit vector + + %% 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 + is_face_hit = false; + hit_position = [NA, NA]; + hit_distance = Inf; + return; + end + + % calculating hit position + hit_position=face.vertex1+kf*tf; + hit_distance = tb; % distance along the light beam +end + diff --git a/beam_trace.m b/beam_trace.m new file mode 100644 index 0000000..401391a --- /dev/null +++ b/beam_trace.m @@ -0,0 +1,31 @@ +function ret = beam_trace(beams, faces, borders) + %% 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 + + 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) + beam_trace(new_beams, faces, borders ); + continue; + end + + % 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 + end +end + + diff --git a/coupling_angles.m b/coupling_angles.m index df7c5f1..ca7dcf0 100644 --- a/coupling_angles.m +++ b/coupling_angles.m @@ -144,4 +144,5 @@ text(0,1.2, str); axis([-1,1,-0.5,1.5],'equal'); %% output of the plot to the file -print('prism_disk_coupling.pdf'); +print('prism_disk_coupling.eps','-depsc2'); +print('prism_disk_coupling.png'); diff --git a/example.m b/example.m new file mode 100644 index 0000000..0b1e855 --- /dev/null +++ b/example.m @@ -0,0 +1,30 @@ +beam.k=[1,-.5]; +beam.origin=[0,0]; +beam.face=NA; + +face.vertex1=[1.1,-1]; +face.vertex2=[2.7,1]; +face.n_right = 1; +face.n_left= 2.5; + +face2.vertex1=[-.1,-3]; +face2.vertex2=[-.1,3]; +face2.n_right = 1; +face2.n_left= 1; + +face3.vertex1=[2.1,-3]; +face3.vertex2=[2.1,3]; +face3.n_right = 1.5; +face3.n_left= 1; + +faces={face,face2,face3}; + +[is_face_hit, hit_position, hit_distance, new_beams] = face_beam_interaction(beam, faces); +beam_reflected=new_beams{1}; +beam_refracted=new_beams{2}; + + +[is_face_hit, hit_position, hit_distance, new_beams] = face_beam_interaction(beam_reflected, faces); + +[is_face_hit, hit_position, hit_distance, new_beams] = face_beam_interaction(beam_refracted, faces); + diff --git a/example_beam_trace.m b/example_beam_trace.m new file mode 100644 index 0000000..ada35c9 --- /dev/null +++ b/example_beam_trace.m @@ -0,0 +1,45 @@ +beam.k=[1,-.5]; +beam.origin=[0,0]; +beam.face=NA; + +face.vertex1=[1.1,-1]; +face.vertex2=[2.7,1]; +face.n_right = 1; +face.n_left= 2.5; + +face2.vertex1=[-.1,-3]; +face2.vertex2=[-.1,3]; +face2.n_right = 1; +face2.n_left= 1; + +face3.vertex1=[2.1,-3]; +face3.vertex2=[2.1,3]; +face3.n_right = 1.5; +face3.n_left= 1; + +faces={face,face2,face3}; + + +border1.vertex1=[-5,-5]; +border1.vertex2=[5,-5]; +border1.n_right = 1.; +border1.n_left= 1; + +border2.vertex1=[5,-5]; +border2.vertex2=[5,5]; +border2.n_right = 1.; +border2.n_left= 1; + +border3.vertex1=[5,5]; +border3.vertex2=[-5,5]; +border3.n_right = 1.; +border3.n_left= 1; + +border4.vertex1=[-5,5]; +border4.vertex2=[-5,-5]; +border4.n_right = 1.; +border4.n_left= 1; + +borders={border1,border2, border3, border4}; +beams={beam}; +beam_trace(beams, faces, borders ) diff --git a/face_beam_interaction.m b/face_beam_interaction.m new file mode 100644 index 0000000..4294c11 --- /dev/null +++ b/face_beam_interaction.m @@ -0,0 +1,107 @@ +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 ibeam starts from face then its index is here + % 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 - index of refraction on the left hand side + % with respect to 1st -> 2nd vertex direction + % face.n_right - index of refraction on the right hand side + + + k=beam.k; + + %% 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 + t=linspace(0,1); + x=face.vertex1(1)+kf(1)*t; + y=face.vertex1(2)+kf(2)*t; + plot(x,y,'k-'); + t=linspace(0,hit_distance); + % draw beam + x=beam.origin(1)+k(1)*t; + y=beam.origin(2)+k(2)*t; + plot(x,y,'r-'); + + % 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; + n2=face.n_right; + else + % beam coming from the right + n1=face.n_right; + n2=face.n_left; + 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); + + % reflected beam direction + theta_normal = atan2(nf(2), nf(1)); + theta_reflected = theta_normal + pi - theta_i; + + beam_reflected.origin = hit_position; + beam_reflected.k = [cos(theta_reflected), sin(theta_reflected)]; + beam_reflected.face=closest_face_index; + new_beams{1} = beam_reflected; + + + % refracted beam direction + % refracted angle with respect to normal + sin_theta_refracted_rel2normal = n1/n2*sin(theta_i); + if ( abs(sin_theta_refracted_rel2normal) >=1 ) + % total internal reflection + else + % beam refracts + theta_refracted_rel2normal = asin( sin_theta_refracted_rel2normal ); + theta_refracted = theta_normal + theta_refracted_rel2normal; + + beam_refracted.origin = hit_position; + beam_refracted.k = [cos(theta_refracted), sin(theta_refracted)]; + beam_refracted.face=closest_face_index; + new_beams{2} = beam_refracted; + end +end + + + -- cgit v1.2.3