% calc_ellipticity
% 121008 TDM
%input is state vector: state = [Jx Jy]

Jx          = state(1);
Jy          = state(2);

Ax          = abs(Jx);
Ay          = abs(Jy);

phi         = angle(Jy)-angle(Jx);

abs_phase   = angle(Jx);

alpha       = atan2(Ay,Ax);
epsilon     = asin(sin(2*alpha)*sin(phi))/2;
theta       = atan(tan(2*alpha)*cos(phi))/2;

disp(['Rotation angle = ' num2str(theta*180/pi) ' degrees '])
disp(['Ellipticity = ' num2str(epsilon*180/pi) ' degrees'])
if epsilon < 0 
    disp('Polarization displays RHC characteristics.')
elseif epsilon > 0 
    disp('Polarization displays LHC characteristics.')
else
    disp('Polarization is linear.')
end