t = linspace(0,2*pi,1000);

phi = pi/2;

E0x = 1;

E0y = 5*sqrt(3);



x = E0x*cos(t);

y = E0y*cos(t+phi);

plot(x,y);



% Principal axes



E0 = sqrt(E0x^2+E0y^2);

ex = 1/E0*[E0x; E0y];

ey = 1/E0*[-E0x; E0y];

xp = E0*cos(phi/2)*cos(t+phi/2);

yp = E0*sin(phi/2)*-sin(t+phi/2);



N = length(t);

XY = zeros(2,N);

for i = 1:N

    XY(:,i) = xp(i)*ex+yp(i)*ey;

end

hold on;

plot(XY(1,:),XY(2,:),'r--');