%% Ce code illustre le transport parallèle d'un vecteur sur une sphère.
%
% (c) Olivier Pujol (Sept. 2014)

close all
clc
clear all

%% Angles des coordonnées sphériques (en radians)
theta=(60:-30:0)*pi/180;    % Colatitude
phi=(0:20:360)*pi/180;      % Azimut
ang1=90;                    % Angle initial entre le vecteur V à transporter et le vecteur unitaire local e_\varphi

%% Initialisation des vecteurs de sortie (colonne 1 = coordonnées cartésiennes X,Y,Z. 
% Colonnes 2 et 3 = angles theta et phi respectivement).

CoordA=nan(3,length(theta),length(phi));    % Point A sur la sphère
V=nan(3,length(theta),length(phi));         % Vecteur V (à transporter)

%% Calcul (transport parallèle)
for iphi=1:length(phi)              % Boucle sur les méridiens (azimut phi fixé)
    ang=(ang1*pi/180)-phi(iphi);    % Modification de l'angle en changeant de méridien
    for itheta=1:length(theta)      % Boucle sur la colatitude
        % Composante du vecteur V dans la base sphérique locale
        V(1,itheta,iphi)=-cos(ang).*sin(phi(iphi))-sin(ang).*cos(phi(iphi)).*cos(theta(itheta));
        V(2,itheta,iphi)=cos(ang).*cos(phi(iphi))-sin(ang).*sin(phi(iphi)).*cos(theta(itheta));
        V(3,itheta,iphi)=sin(ang).*sin(theta(itheta));
        % Coordonnées sphériques du point A
        CoordA(1,itheta,iphi)=cos(phi(iphi)).*sin(theta(itheta));
        CoordA(2,itheta,iphi)=sin(phi(iphi)).*sin(theta(itheta));
        CoordA(3,itheta,iphi)=cos(theta(itheta));
    end
end
     
V=0.5*V;    % Réduction de la longueur du vecteur

%% Tracé

% Première présentation (figure 1)

figure(1)
[x,y,z]=sphere(18);
h=surf(x,y,z);
alpha(0.25);
axis off;axis equal
colormap(bone)
shading faceted
set(h,'LineStyle',':','LineWidth',0.5)
hold on

colat=0:pi/50:theta(1); % Vecteur des colatitudes (pour tracé ultérieur des méridiens) 

% Transport parallèle de H1 à N (phi = 0 ; iphi = 1)
% Tracé des cercles de colatitude
az=-pi/2+phi(1)+20*pi/180:pi/50:pi/2+phi(1)+20*pi/180;
h1=plot3(cos(az),sin(az),0.*az,'k','linewidth',1);
h2=plot3(sin(theta(1))*cos(az),sin(theta(1))*sin(az),cos(theta(1)).*ones(size(az)),'k','linewidth',1);
h3=plot3(-cos(az),-sin(az),0.*az,'k--','linewidth',1);
h4=plot3(-sin(theta(1))*cos(az),-sin(theta(1))*sin(az),cos(theta(1)).*ones(size(az)),'k--','linewidth',1);

for itheta=1:length(theta)
    figure(1)
    title(['Azimut $\varphi=',num2str(180*phi(1)/pi), '^\circ$, '...
        'colatitude $\theta=',num2str(180*theta(itheta)/pi), '^\circ$'],...
        'Interpreter','latex','Fontsize',12)
    view(110+phi(1)*180/pi,25);
    % Tracé des méridiens
    plot3(sin(colat)*cos(phi(1)),sin(colat)*sin(phi(1)),cos(colat),'k','linewidth',1);
    % Tracé du vecteur transporté
    quiver3(CoordA(1,itheta,1),CoordA(2,itheta,1),CoordA(3,itheta,1),...
          V(1,itheta,1),V(2,itheta,1),V(3,itheta,1),...
            'b','linewidth',2,'Marker','.','Markersize',10,'MarkerEdgeColor','k');
    pause(0.5)
    
    %Supprimer le commentaires de la ligne suivante pour imprimer les figures en eps
    %print(figure(1),'-depsc2',['TraPar_az',num2str(iphi),'.eps']); %
    
end

% Transport parallèle de N à H2 (pas = 60°)
for iphi=4:3:length(phi) 
    delete(h1,h2,h3,h4)
    clear az h1 h2 h3 h4
    % Tracé des cercles de colatitude
    az=-pi/2+phi(iphi)+20*pi/180:pi/50:pi/2+phi(iphi)+20*pi/180;
    h1=plot3(cos(az),sin(az),0.*az,'k','linewidth',1);
    h2=plot3(sin(theta(1))*cos(az),sin(theta(1))*sin(az),cos(theta(1)).*ones(size(az)),'k','linewidth',1);
    h3=plot3(-cos(az),-sin(az),0.*az,'k--','linewidth',1);
    h4=plot3(-sin(theta(1))*cos(az),-sin(theta(1))*sin(az),cos(theta(1)).*ones(size(az)),'k--','linewidth',1);

    for itheta=length(theta):-1:1
       figure(1)
       title(['Azimut $\varphi=',num2str(180*phi(iphi)/pi), '^\circ$, '...
        'colatitude $\theta=',num2str(180*theta(itheta)/pi), '^\circ$'],...
        'Interpreter','latex','Fontsize',12)
       view(110+phi(iphi)*180/pi,25);
       % Tracé des méridiens
       plot3(sin(colat)*cos(phi(iphi)),sin(colat)*sin(phi(iphi)),cos(colat),'k','linewidth',1);
       quiver3(CoordA(1,itheta,iphi),CoordA(2,itheta,iphi),CoordA(3,itheta,iphi),...
            V(1,itheta,iphi),V(2,itheta,iphi),V(3,itheta,iphi),...
            'b','linewidth',1.5,'Marker','+','Markersize',5,'MarkerEdgeColor','k');hold on
       pause(0.5)      
    end
    
    % Vecteur initial (phi=0)
    Vinit=[cos(phi(iphi)) sin(phi(iphi)) 0;sin(phi(iphi)) cos(phi(iphi)) 0;0 0 1]*V(:,1,1);
    quiver3(CoordA(1,1,iphi),CoordA(2,1,iphi),CoordA(3,1,iphi),...
          Vinit(1,1,1),Vinit(2,1,1),Vinit(3,1,1),'k--','linewidth',1.5);
    pause(0.5)
       
    %Supprimer le commentaires de la ligne suivante pour imprimer les figures en eps
    %print(figure(1),'-depsc2',['TraPar_az',num2str(iphi),'.eps']);
end

% Autre représentation (figure 2)

figure(2)
[x,y,z]=sphere(18);
h=surf(x,y,z);
view(90+phi(1)*180/pi,90-theta(1)*180/pi);
alpha(0.25);
axis equal;axis off
colormap(bone)
shading faceted
set(h,'LineStyle',':','LineWidth',0.5)
hold on

az=-pi/2:pi/50:pi/2;
plot3(cos(az),sin(az),0.*az,'k','linewidth',1);
plot3(-cos(az),sin(az),0.*az,'k--','linewidth',1);
plot3(sin(theta(1))*cos(az),sin(theta(1))*sin(az),cos(theta(1)).*ones(size(az)),'k','linewidth',1);
plot3(-sin(theta(1))*cos(az),sin(theta(1))*sin(az),cos(theta(1)).*ones(size(az)),'k--','linewidth',1);

for iphi=1:length(phi)
    figure(2)
    hold all
    hq=quiver3(CoordA(1,1,iphi),CoordA(2,1,iphi),CoordA(3,1,iphi),...
        V(1,1,iphi),V(2,1,iphi),V(3,1,iphi),1,...
        'b','linewidth',1.5,'Marker','+','Markersize',5,'MarkerEdgeColor','k');
    pause(0.5)
end

%Supprimer le commentaires de la ligne suivante pour imprimer les figures en eps
%print(figure(2),'-depsc2','TraPar.eps');
