% Returns the reflectivities as a function of incident angle. rs and rp are
% NxM matrixes, each row corresponds to a particular wavelength and each column
% to a particular angle (i.e. wavelength is labeled by N and angle by M).
%
% na = Nx1, Nx2 or Nx3 matrix of indexes of the medium on the near side of the interface (i.e. the ambient medium).
%      May be given as a 1xm vector (m=1,2,3) if nb does not change with na.
%      na(N,1), na(N,2) and na(N,3) are the isotropic indices as defined in Orfanidi's function fresnel.m
%      (which is called by this function). The firstindex would often label the indices at different wavelengths.% nb = Nx1, Nx2 or Nx3 matrix of indexes of the medium on the far side of the interface.
%      May be given as a 1xm vector (m=1,2,3) if na does not change with nb.
%      nb(N,1), nb(N,2) and nb(N,3) are the isotropic indices as defined in Orfanidi's function fresnel.m
%      (which is called by this function). The firstindex would often label the indices at different wavelengths.
% theta = vector of angles of incidence in degrees.
%
% Relies on the function fresnel.m by Orfanidis which must be in the path.
% If only 1-D matrixes are needed and the medium is isotropic, use fresneliso.m instead.  It is much faster.

function [rs rp]=fresnel2(na,nb,theta)

if size(na,1)==1
    na=repmat(na,size(nb));  
elseif size(nb,1)==1
    nb=repmat(nb,size(na));
end

rs=zeros(size(nb,1),length(theta)); rp=rs;

for s=1:size(nb,1)
    [rs(s,:) rp(s,:)]=fresnel(na(s,:),nb(s,:),theta);
end