function [khat_r, xhat_p, yhat_pi, yhat_pr, MP1, MP2, a1, J2, a2] ...
    = find_global_reflection_from_surface(J1,N_i,N_t,khat_i,nhat);

%Note:  J1 is assumed to be defined with xhat_p and yhat_pi. 

%find reflected wave direction
khat_r      = -2*dot(khat_i,nhat)*nhat + khat_i;

%find orientation of local jones axes
if norm(cross(khat_r,khat_i))~=0
    %calculate incident xhat_p and yhat_p
    xhat_p      = cross(khat_r,khat_i);
    xhat_p      = xhat_p/norm(xhat_p);
else  %k vectors are either parallel or antiparallel
    xhat_p      = cross(khat_i,circshift(khat_i,[0 1]));
    xhat_p      = xhat_p/norm(xhat_p);
end
%calculate incident yhat_p
yhat_pi     = cross(khat_i,xhat_p);
%calculate reflected yhat_p
yhat_pr     = cross(khat_r,xhat_p);

%find s and p components of incident light
%are consistent with Fresnel reflection convention
%introduced in Ch 3 (2009).
shat            = xhat_p;   %we don't really use these, 
                            %but they are listed for clarity
phat_i          = yhat_pi;
phat_r          = yhat_pr;

%Calculate cosines
cos_theta_i = abs(dot(khat_i,nhat));
cos_theta_t = (1-(N_i/N_t)^2 * (1-cos_theta_i.^2))^(0.5);

%Calculate reflectivities
r_s         = (N_i*cos_theta_i-N_t*cos_theta_t)/(N_i*cos_theta_i+N_t*cos_theta_t);
r_p         = (N_t*cos_theta_i-N_i*cos_theta_t)/(N_t*cos_theta_i+N_i*cos_theta_t);

%Apply mirror matrix
M_m         = [r_s 0;0 r_p];
J2          = M_m*J1;

%find MP1 and MP2
MP1         = [xhat_p(:) yhat_pi(:)];
MP2         = [xhat_p(:) yhat_pr(:)];

%find a1 and a2:
a1          = MP1*J1;
a1          = a1/norm(a1);
a2          = MP2*J2;
a2          = a2/norm(a2);