% Dominic Grenier octobre 2017
% Trace l'orientation du champ E en fonction du temps sur T/4 en rouge
% direction de propagation a_P (ou a_beta) en bleu
% sens de rotation de la polarisation en vert
% 
% En utilisant "Rotate 3D" de la fenêtre graphique, on peut:
% se placer dans le plan E pour voir si E est bien transverse à a_P
% s'aligner avec l'axe de propagation a_P pour voir la polarisation
%
% Rotation des axes pour aligner l'axe de propagation avec z' en cyan
% de sorte que le champ E en magenta se trouve dans le plan x'-y'
% sens de rotation de la propagation en jaune dans système x'y'z'

% Exemple avec Q.17 Chap.5
% E0=[3+j*0; 4+j*0; 0+j*5];
% a_beta=[0.8; -0.6; 0];

E=[3+j*0; 4+j*5; 0+j*10];
a_beta=[0.8; -0.6; 0.3];


% E(x,y,z,t=0:0.05:0.25)
% début du code
figure(1)
mm=max(abs(E));
beta=mm;
omega=2*pi*1;
%T=[0,0.05,0.25]';
T=0:0.05:0.5; lT=length(T);

Ei=(abs(E)*ones(1,lT)).*cos(omega*(ones(3,1)*T)+(angle(E)*ones(1,lT)));

ma_beta=sqrt(a_beta'*a_beta);
Betai=[mm*a_beta/ma_beta];

theta=sign(a_beta(3))*acos(a_beta(3)/ma_beta);
phi=atan2(a_beta(2),a_beta(1));
Et=[cos(phi) sin(phi) 0; -sin(phi) cos(phi) 0; 0 0 1]*Ei;
Ef=[cos(theta) 0 -sin(theta);0 1 0; sin(theta) 0 cos(theta)]*Et;
Betat=[cos(phi) sin(phi) 0; -sin(phi) cos(phi) 0; 0 0 1]*Betai;
Betaf=[cos(theta) 0 -sin(theta);0 1 0; sin(theta) 0 cos(theta)]*Betat;

Ex0=max(Ef(1,:)); Ey0=max(Ef(2,:));
alpha= atan(Ey0/Ex0)/d2r;
Emax=sqrt(max(Ef(1,:).^2+Ef(2,:).^2));
Emin=sqrt(min(Ef(1,:).^2+Ef(2,:).^2));
AR=Emax/Emin; fprintf('rapport=%6.2f\n',AR);

Or=zeros(length(T),1);
quiver3(Or,Or,Or,Ei(1,:)',Ei(2,:)',Ei(3,:)',1,'r');
hold on;
xlabel('x');ylabel('y');zlabel('z');
axis([-mm,mm,-mm,mm,-mm,mm]); axis equal;

dEx=gradient(Ei(1,:));
dEy=gradient(Ei(2,:));
dEz=gradient(Ei(3,:));
quiver3(Ei(1,1:end-1),Ei(2,1:end-1),Ei(3,1:end-1),dEx(1:end-1),dEy(1:end-1),dEz(1:end-1),'g');
hold on
quiver3(0,0,0,Betai(1),Betai(2),Betai(3),'b','LineWidth',2,'MaxHeadSize',0.25);
hold on

quiver3(Or,Or,Or,Ef(1,:)',Ef(2,:)',Ef(3,:)',1,'m');
dEx=gradient(Ef(1,:));
dEy=gradient(Ef(2,:));
Or=zeros(1,length(T)-1);
quiver3(Ef(1,1:end-1),Ef(2,1:end-1),Or,dEx(1:end-1),dEy(1:end-1),Or,'y');
hold on
quiver3(0,0,0,Betaf(1),Betaf(2),Betaf(3),'c','LineWidth',2,'MaxHeadSize',0.25);
hold on

axis equal
hold off;