%% This script aims at illustrate vectorial parallel transport 
%% over a sphere.
%
% (c) Olivier Pujol (Sept. 2014)

close all
clc
clear all

%% Spherical angles (in radians)
theta=(60:-30:0)*pi/180;    % Colatitude
phi=(0:20:360)*pi/180;      % Azimuth
ang1=90;                    % Angle between initial vector V and the local varphi-eigenvector

%% Output vectors initialisation (first column is cartesian coordinates X,Y,Z. 
% The second and third ones are respectively the values of theta and phi

CoordA=nan(3,length(theta),length(phi));    % Point A over a sphere
V=nan(3,length(theta),length(phi));         % Vector V (to be transported)

%% Calculation (parallel transport)
for iphi=1:length(phi)              % Loop over the meridian (azimuth phi)
    ang=(ang1*pi/180)-phi(iphi);    % Angle modification when changing meridian
    for itheta=1:length(theta)      % Loop over the colatitude
        % Component of vector V in the local spherical basis
        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));
        % Spherical coordinate of 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;    % Factor of reduction to limit vector length

%% Plotting

% First presentation (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 of colatitudes (for further meridian lines tracing) 

% Parallel transport from H1 to N (phi = 0 ; iphi = 1)
% Latitude lines tracing
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(['Azimuth $\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);
    % Meridian line tracing
    plot3(sin(colat)*cos(phi(1)),sin(colat)*sin(phi(1)),cos(colat),'k','linewidth',1);
    % Trace of the displaced vector
    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)
    
    %Uncomment the following line to print figure in eps
    %print(figure(1),'-depsc2',['TraPar_az',num2str(iphi),'.eps']); %
    
end

% Parallel transport from N to H2 (step = 60°)
for iphi=4:4%3:length(phi) 
    delete(h1,h2,h3,h4)
    clear az h1 h2 h3 h4
    % Latitude lines tracing
    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(['Azimuth $\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);
       % Meridian lines tracing
       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
    
    % Initial vector (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)
       
    %Uncomment the following line to print figure 1 in eps
    %print(figure(1),'-depsc2',['TraPar_az',num2str(iphi),'.eps']);
end
pause
% Other representation (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

%Uncomment the following line to print figure 2 in eps
%print(figure(2),'-depsc2','TraPar.eps');
