beam1.k=[-1,0];
beam1.origin=[2,0.5];
%beam1.k=[0,-1];
%beam1.origin=[.4,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.4885, 2.7532]; % n_s and n_p

face1.vertex1=[-1,0];
face1.vertex2=[1,0];
face1.n_right = [1, 1];
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');