function integration_at_L_prime
%integration uses s only, no derivative wrt lambda etc

clear; clf;
L=10;                                        %max x range
L_prime=-L;
lambda_0=-2.0-0.5*i;                               %centre of lambda path circle
r=0.8;                                        %radius of lambda path circle

%setup graphics boxes
eigenvalues=[-0.75,0,1.25];
subplot(131); xlabel Re(\lambda); ylabel Im(\lambda);
grid on; hold on; axis([-3.0,5,-4,2]);
plot(real(eigenvalues), imag(eigenvalues),'*k');
title(['\lambda -path: \lambda_0 = ', num2str(lambda_0), ' , r = ',num2str(r)]);
xlabel Re(\lambda(s)); ylabel Im(\lambda(s)); hold on; grid on;

numsteps = 100; %(just for the picture)
s= 0:2*pi/numsteps:2*pi;
lambda_array = lambda_0 + r*exp(i*s);
subplot(131)
plot(real(lambda_array), imag(lambda_array),'r');

subplot(132); title('s-\theta/(2\pi) plot'); xlabel s; ylabel \theta(s)/(2\pi); hold on; grid on;
subplot(131);
f=getframe;

%integrate the connection around gamma_3
theta_init=0;
j=1;
s_numsteps=10000;
delta = 2*pi/s_numsteps;
sspan=0:0.05:2*pi;%[0 2*pi];
[s theta_array] = ode45(@spiral_theta_integration, sspan, theta_init);
subplot(132)
theta_array=theta_array/(2*pi);
plot(s,theta_array)
delta_theta = theta_array(length(theta_array));
title(['s-\theta/(2\pi) plot, \Delta(\theta)/2\pi=',num2str(delta_theta)])


%S^2 figure;
subplot(133)
[a b]=meshgrid(-1:0.05:1,-1:0.05:1);z=real(sqrt(1-a.^2-b.^2));
surf(a,b,z); hold on; title('Image of path on S^3 under Hopf map');
z=-real(sqrt(1-a.^2-b.^2));surf(a,b,z);axis equal;set(gca,'ZDir','reverse');
xlabel '\pi_1(v)'; ylabel '\pi_2(v)'; zlabel '\pi_3(v)';
plot3(path_on_S2_array(:,1),path_on_S2_array(:,2),path_on_S2_array(:,3));
comet3(path_on_S2_array(:,1),path_on_S2_array(:,2),path_on_S2_array(:,3));


a=length(theta_array); %new a,b
b=length(path_on_S2_array);
if a>=b
    for j=1:a
        k=ceil((j/a)*b);
        subplot(131); plot(real(lambda_0 + r*exp(i*s(j))),imag(lambda_0 + r*exp(i*s(j))),'.k');
        subplot(132); plot(s(j),theta_array(j),'.k');
        subplot(133); plot3(path_on_S2_array(k,1),path_on_S2_array(k,2),path_on_S2_array(k,3),'.y');
        F(j)=getframe(gcf);
    end
else
    for j=1:b
        k=ceil((j/b)*a);
        subplot(131); plot(real(lambda_0 + r*exp(i*s(k))),imag(lambda_0 + r*exp(i*s(k))),'.k');
        subplot(132); plot(s(k),theta_array(k),'.k');
        subplot(133); plot3(path_on_S2_array(j,1),path_on_S2_array(j,2),path_on_S2_array(j,3),'.y');
        F(j)=getframe(gcf);
    end
end

%pause;
%figure;
%movie(F);
%pause;
%movie2avi(F,'integrate_at_Lprime_0.avi');

%Nested Functions
%----------------------------------------------------------
%function for integrating along spiral path
    function theta_s = spiral_theta_integration(s,theta)
        
        present_s = s;
        next_s = s + delta; 
        last_s = s - delta;
        
        present_lambda = lambda_0 + r*exp(i*present_s);   
        next_lambda=lambda_0 + r*exp(i*(next_s));          
        last_lambda=lambda_0 + r*exp(i*(last_s));
        
        present_mu=sqrt(present_lambda+1);
        next_mu=sqrt(next_lambda+1);
        last_mu=sqrt(last_lambda+1);

        %next_lambda... integration of u from x=L -> L_prime
        xspan=[L L_prime];           %x integration range
        lambda=next_lambda; mu=next_mu;
        ksi_1=[ 1; -mu]; u_init=ksi_1;
        [x u] = ode45(@u_ODE, xspan, u_init);
        next_u=u(length(u),:).'; %results of integration
        clear u

        %present_lambda... integration of u from x=L -> L_prime
        xspan=[L L_prime];
        lambda=present_lambda; mu=present_mu;
        ksi_1=[ 1; -mu]; u_init=ksi_1;
        [x u] = ode45(@u_ODE, xspan, u_init);
        present_u=u(length(u),:).';
        clear u
       
        %last_lambda... integration of u from x=L -> L_prime
        xspan=[L L_prime];
        lambda=last_lambda; mu=last_mu;
        ksi_1=[ 1; -mu]; u_init=ksi_1;
        [x u] = ode45(@u_ODE, xspan, u_init);
        last_u=u(length(u),:).';
        clear u
       
        u=present_u;
        
        mod_u=sqrt(u'*u);
        v=u/mod_u;                                               %v=(v1,v2) in S^3 in C^2
       
        w(3)=conj(v(1))*v(2) + v(1)*conj(v(2));
        w(2)=i*(conj(v(1))*v(2) - v(1)*conj(v(2)));
        w(1)= (abs(v(1)))^2 - (abs(v(2)))^2;                        %w=(w1,w2,w3) in S^2 in R^3
        
        path_on_S2_array(j,:)=w;
        j=j+1;
        
        u_s = (next_u - last_u)/(2*delta);

        theta_s = imag((u')*(u_s))/(u'*u);

        %----------------------------------------------------------
        function u_x = u_ODE(x,u)
            f = (3/2)*((sech(x/2))^2);
            A = [	0               1
                    lambda+1-2*f    0	];
            %u_x = A*u;
            u_x = (A - (-mu)*eye(2))*u;
        end
        %----------------------------------------------------------
    end

%----------------------------------------------------------    

end