function [Ex,Ey,psid,AR]=polaris(Ex0,Ey0,dephase,omega,t,beta,z)
% Dominic Grenier octobre 2014

d2r=pi/180;


if (nargin==3)
    Ex=Ex0; Ey=Ey0;
    alpha= atan(Ey0/Ex0);
%    psi = 0.5*atan(tan(2*alpha)*cos(dephase*d2r))
    psi = 0.5*atan2(2*Ex0.*Ey0*cos(dephase*d2r),Ex0.^2-Ey0.^2) %mieux pour considérer E_y>E_x

    xi = 0.5*asin(2*Ex0*Ey0*sin(dephase*d2r)/(Ex0^2+Ey0^2))
%    xi  = 0.5*asin(sin(2*alpha)*sin(dephase*d2r))
    
    AR = abs(1/tan(xi))
%    AR   = abs(Ey0*sin(psi)*sin(dephase*d2r))/abs(Ex0*sin(psi)-Ey0*cos(psi)*cos(dephase*d2r))
%    AR  = abs(1./tan(0.5*asin(sin(2*alpha)*sin(dephase*d2r))))
    if (AR<1)
        AR = 1/AR;
    end
    psid  = psi/d2r;
else
    alpha= atan(Ey0/Ex0);
%    psi  = 0.5*atan(tan(2*alpha)*cos(dephase*d2r));
    psi = 0.5*atan2(2*Ex0.*Ey0*cos(dephase*d2r),Ex0.^2-Ey0.^2);
    if psi<0 psi=psi+pi; end
    psid  = psi/d2r;

    xi = 0.5*asin(2*Ex0*Ey0*sin(dephase*d2r)/(Ex0^2+Ey0^2));
    xid = xi/d2r;
    AR = abs(1/tan(xi));
    
    Ex = Ex0*cos(omega*t-beta*z);
    Ey = Ey0*cos(omega*t-beta*z+dephase*d2r);
end

