function [Rs, Rp] = Fresnel(theta,n1,n2)

rooted_term=sqrt(1-(n1*sin(theta)/n2).^2);
R_s = (n1*cos(theta)-n2*rooted_term)./(n1*cos(theta)+n2*rooted_term);
R_p = (n1*rooted_term-n2*cos(theta))./(n1*rooted_term+n2*cos(theta));
Rs = R_s.*conj(R_s);
Rp = R_p.*conj(R_p);

return