function integration_at_L_prime_HSpulse
%integrates connection form around loop of solns at L_prime

clear; clf;
L=10;                                        %max x range
L_prime=-L;
lambda_0=0+0.1*i%;0.7*i;                               %centre of lambda path circle
r=1;                                        %radius of lambda path circle

%HS pulse params
omega = 0.03;%2+(sqrt(5));
rho = 1;%1/(sqrt(5));
psi = pi/2;%atan(2);


%setup graphics boxes
subplot(121); xlabel Re(\lambda); ylabel Im(\lambda);
grid on; hold on; %axis([-10,10,-5,5]);
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(121)
plot(real(lambda_array), imag(lambda_array),'r');

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

%integrate the connection around gamma_3
theta_init=0;
s_numsteps=10000;
delta = 2*pi/s_numsteps;
sspan=0:0.05:2*pi;%[0 2*pi];
[s theta_array] = ode45(@theta_integration, sspan, theta_init);
subplot(122)
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)])


a=length(theta_array); %new a,b
for j=1:a
    subplot(121); plot(real(lambda_0 + r*exp(i*s(j))),imag(lambda_0 + r*exp(i*s(j))),'.k');
    subplot(122); plot(s(j),theta_array(j),'.k');
    F(j)=getframe(gcf);
end

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

%Nested Functions
%----------------------------------------------------------
%function for integrating connection form around solns at L_prime
    function theta_s = theta_integration(s,theta)
        s
        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));

        lambda = present_lambda;
        p = 2*omega + lambda*rho*(sin(psi));
        tau = 1 - omega^2 + lambda*rho*cos(psi);
        sigma_plus = -(sqrt(2))*(sqrt(tau + sqrt(tau^2 + p^2))); %eval of A_infinty of most -ve real part
        present_sigma_plus = sigma_plus;
        present_ksi_plus = [    2*sigma_plus
                                -2*p
                                sigma_plus^2
                                -(sigma_plus^2)
                                -2*p
                                sigma_plus*(sigma_plus^2 - 2*tau)   ];  %evec of A_infinty for sigma_plus
        
        lambda = next_lambda;
        p = 2*omega + lambda*rho*(sin(psi));
        tau = 1 - omega^2 + lambda*rho*cos(psi);
        sigma_plus = -(sqrt(2))*(sqrt(tau + sqrt(tau^2 + p^2))); %eval of A_infinty of most -ve real part
        next_sigma_plus = sigma_plus;
        next_ksi_plus = [   2*sigma_plus
                            -2*p
                            sigma_plus^2
                            -(sigma_plus^2)
                            -2*p
                            sigma_plus*(sigma_plus^2 - 2*tau)   ];  %evec of A_infinty for sigma_plus
        
        lambda = last_lambda;
        p = 2*omega + lambda*rho*(sin(psi));
        tau = 1 - omega^2 + lambda*rho*cos(psi);
        sigma_plus = -(sqrt(2))*(sqrt(tau + sqrt(tau^2 + p^2))); %eval of A_infinty of most -ve real part
        last_sigma_plus = sigma_plus;
        last_ksi_plus = [   2*sigma_plus
                            -2*p
                            sigma_plus^2
                            -(sigma_plus^2)
                            -2*p
                            sigma_plus*(sigma_plus^2 - 2*tau)   ];  %evec of A_infinty for sigma_plus

        %next_lambda... integration of u from x=L -> L_prime
        xspan = [L L_prime];           %x integration range
        lambda = next_lambda;
        sigma_plus = next_sigma_plus;
        u_init = next_ksi_plus;
        [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;
        sigma_plus = present_sigma_plus;
        u_init = present_ksi_plus;
        [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;
        sigma_plus = last_sigma_plus;
        u_init = last_ksi_plus;
        [x u] = ode45(@u_ODE, xspan, u_init);
        last_u=u(length(u),:).';
        clear u
       
        u=present_u;
        
        u_s = (next_u - last_u)/(2*delta);

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

        %----------------------------------------------------------
        function u_x = u_ODE(x,u)

        q1 = sech(x)*cos(omega*log(cosh(x)));
        q2 = -sech(x)*sin(omega*log(cosh(x)));
        a_31 = lambda*rho*cos(psi) + 1 - omega^2 - (2 - omega^2)*(q2^2 + 3*q1^2) + 6*omega*q1*q2;
        a_32 = -lambda*rho*sin(psi) - 2*omega - 2*(2 - omega^2)*q1*q2 + 3*omega*(q1^2 + 3*q2^2);
        a_41 = lambda*rho*sin(psi) + 2*omega - 2*(2 - omega^2)*q1*q2 - 3*omega*(3*q1^2 +q2^2);
        a_42 = lambda*rho*cos(psi) + 1 - omega^2 - (2 - omega^2)*(q1^2 + 3*q2^2) - 6*omega*q1*q2;

        A = [	0		0		1		-1		0		0
                a_32	0		0		0		0		0
                a_42	0		0		0		0		1       %A(x,lambda), which varies along the integration path
               -a_31	0		0		0		0      -1
               -a_41	0		0		0		0		0
                0      -a_41	a_31   -a_42	a_32	0	];

           %u_x = A*u;
            u_x = (A - (sigma_plus)*eye(6))*u;
        end
        %----------------------------------------------------------
    end

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

end