Determining relative angular orientation of vectors

I am trying to implement Snell's law, but I am currently lacking an elegant way to take into account the angular orientation/direction into which I need to rotate my vector. Quick and dirty example:

%choose incident vector and intersection points with circle
vector = [1 ; 0];
P1 = [-sqrt(0.0075); 0.05];
P2 = [-sqrt(0.0075); -0.05];

%for both points, calculate the angle between the incident vector and the
%normal
theta1 = real(acos(max(min(dot(vector, -P1) / (vecnorm(vector) * vecnorm(P1))))));
theta2 = real(acos(max(min(dot(vector, -P2) / (vecnorm(vector) * vecnorm(P2))))));

%snell's law
alpha1 = asin((sin(theta1)/1.6));
alpha2 = asin((sin(theta2)/1.6));

%rotate the normal vector by the refracted angle (*10 so the length is the
%same as for version 2)
new_vec1 = rot_matrix_2d(-P1, -alpha1)*10;
new_vec2 = rot_matrix_2d(-P2, alpha2)*10;

%version 2: rotate incident vector by difference between incident and
%refracted angle
new_vec1_v2 = rot_matrix_2d(vector, (theta1-alpha1));
new_vec2_v2 = rot_matrix_2d(vector, -(theta2-alpha2));

%definition of the rotation matrix
function rot_vec = rot_matrix_2d(vec, theta)
    R = [cos(theta) sin(theta) ; -sin(theta) cos(theta)];
    rot_vec = R * vec;
end

I have a vector (or rather two) that goes to the right and want to make it refract on a circular surface at two points (x0 = y0 = 0, r = 0.1). The upper vector should, after refraction, go to the bottom right, while the lower vector should go to the upper right. P1 and P2 are the intersection points of the vector(s) and the circle, and as such, -P1 and -P2 are also the normals. The issue basically starts when determining the angle, since theta1 = theta2 in value. However, in one case, the angle describes a rotation in clockwise direction, and in the other case, it is counter clockwise.

Getting the refracted vector can be done by either rotating the normal by the angle after refraction (here with index of refraction n = 1.6), or by rotating the original vector by the difference between the refracted angle and the angle of incidence (delta). However, as you can see, in one case I need to rotate by -alpha, while in the other case, I need to rotate by alpha (or delta and -delta). Is there an easy way to determine into which direction I need to rotate? I want to make this generically applicable.

[Jan]

Some comments in the code would clarify the purpose.

[Dominik Rhiem]

@Jan added some comments.

[Matt J]

Is there an easy way to determine into which direction I need to rotate?

It's been a long time since physics class. Does Snell's law say which way the line should bend? WIll it always be toward the normal or will it be away from the normal? Which side of the normal will it be on?

[Dominik Rhiem]

@Matt J that depends on the indexes of refraction. If the ray enters a medium with a higher refractive index, it will bend towards the normal, if it enters a medium with a lower refraction index, it will bend away. Here is an example for entering a medium with a higher refractive index.

Note that the ray paths are also correct in the "other" direction, i.e. if you turn the direction of the ray around, this image is still correct.

[Matt J]

Yes, but knowing the angle with the normal isn't enough to resolve the ambiguity. There are two different vectors making the same angle with the normal, as illustrated below:

[Dominik Rhiem]

@Matt J Ah sorry, I misunderstood the question then. To put it this way: the ray "crosses" the normal, it continues on the "other side".

[Matt J]

Then see my answer below.

[Matt J]

To put it this way: the ray "crosses" the normal, it continues on the "other side".

Apparently it is possible to have a negative index of refraction, in which case the beam will not cross the normal:

https://en.wikipedia.org/wiki/Negative-index_metamaterial

Answer 1

[Matt J]

Perhaps as follows?

vector = [1 ; 0];
P1 = [-sqrt(0.0075); 0.05];
P2 = [-sqrt(0.0075); -0.05];

new_vec1 = refract(vector,P1,1.6)

new_vec1 = 2×1

0.9789 -0.2043

new_vec2 = refract(vector,P2,1.6)

new_vec2 = 2×1

0.9789 0.2043

function refractedVector = refract(vector,normal,snellFactor)

   vector=vector(:)/norm(vector);
   normal=normal(:)/norm(normal);


   normal = sign( dot(vector,normal) )*normal; %direct the normal so that incident angle is acute
   
   alpha= asin( sin( acos( dot(vector,normal) )  )   /snellFactor   ); %Snell's law
      

   s=sign( det([vector,normal])); %s<0 ==> clockwise
   
   antinormal=normal([2,1]).*[1;-1].*s;
   
   refractedVector=normal+tan(alpha)*antinormal;
   refractedVector=refractedVector/norm(refractedVector);
   
end

[Dominik Rhiem]

Also here sorry for my late answer. I had to change a few things to vectorize the code, but it works, so thanks!

Answer 2

[Bruno Luong]

No need to call costly trigonometric functions

% Direction of incidence Ray (normalized)
Rayi = randn(2,1);
Rayi = Rayi/norm(Rayi);

% Direction of normal vector  (normalized)
Normal = randn(2,1);
Normal = Normal/norm(Normal);

% Refractive index, n1 on incidence side, n2 on refraction side
n1 = 1;
n2 = 1.6;

Tangent = [0 -1; 1 0]*Normal;
Ci = dot(Rayi,Normal);
Si = dot(Rayi,Tangent);
Sr = (n1/n2)*Si;
if abs(Sr) > 1
    % Should never get here if n2 > n1
    fprintf('Total reflexion: No refraction ray\n')
    return
end
Cr = sign(Ci)*sqrt(1-Sr^2); % NOTE: The sign(Ci) can be removed if N is by convention pointed from n1 to n2 medium.
Rayo = Cr*Normal + Sr*Tangent;

Rayi

Rayi = 2×1

-0.8475 -0.5308

Normal

Normal = 2×1

-0.9958 -0.0916

Rayo

Rayo = 2×1

-0.9296 -0.3685

[Dominik Rhiem]

Sorry, I'm just now getting around to testing out your solution. I like your solution for avoiding the trigonometric functions. How exactly can I understand Sr? You use it to take care of total reflexion, which is very useful since I want to implement reflexion losses, as well. I suppose I can include these by multiplying my amplitude by , is that correct?

[Bruno Luong]

The reflexion direction is (not tested)

Rayr = -Cr*Normal + Si*Tangent;

If you want to compute amplitude of the light refraction/relexion you need to use https://en.wikipedia.org/wiki/Fresnel_equations

[Dominik Rhiem]

@Bruno Luong I had them actually already implemented and still forgot... no way around getting the angles then. Thanks anyway!

[Bruno Luong]

Why you want to insist on getting the angle?

[Dominik Rhiem]

@Bruno Luong because I want to take amplitude reduction due to reflection losses into account, which depends on the involved angles, as shown in the link you posted.

[Bruno Luong]

For amplitude as well as direction, you don't need the angles, all you need are (square of) sine and cosine of different rays, which is given by (incident, transmitted and reflexion):

• cosine(thetai) = Ci
• sine(thetai) = Si
• cosine(thetat) = Cr
• sine(thetat) = Sr
• cosine(thetar) = Ci
• sine(thetat) = Si