function [khi,gamma]=polarization(ax, ay, delta)
% returns the rotation angle, gamma, and the eplipticity angle, khi in degrees
% ax and ay are scalars while delta is the phase angle difference in degrees

deltar=delta*pi/180;

psir=atan(ay/ax);
gammar=1/2*(atan((tan(2*psir)*cos(deltar))));
khir=1/2*asin(sin(2*psir)*sin(deltar));
khi=180*khir/pi;
%convert from radians to degrees
gamma=180*gammar/pi;
if cos(deltar)>0 & gamma<0
   gamma=gamma+90;
elseif cos(deltar)<0 & gamma>0
   gamma=gamma-90;
else
   gamma=gamma;
end



