Fresnel equation vectorized with respect to incident angle

1 files changed, 11 insertions, 10 deletions

diff --git a/fresnel_reflection.m b/fresnel_reflection.m
index 0e6dfdc..347d60b 100644
--- a/fresnel_reflection.m
+++ b/fresnel_reflection.m
@@ -1,27 +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
+%% 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) );
 
-	if ( abs(sin_theta_t) >= 1)
-		%% total internal reflection
-		R = [1, 1];
-		theta_t=sign(theta_i)*pi/2;
-		return;
-	end
 
 	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];
+	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