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

clear; clf;
L=12;%13.5;                                        %max x range
L_prime=-L;
lambda_0=-2;                               %centre of lambda path circle
r=0.4;                                        %radius of lambda path circle

%Schrodinger params
D = 12;
omega = 0.204124;
n = 10:1:23;
E_n = -12 + (n+0.5) - (1/48)*(n+0.5).^2;


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


%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;
        if real(-sqrt(-2*lambda))<real(sqrt(-2*lambda))
        present_ksi_plus=[ -sqrt(-2*lambda); 0];                                  %stable evec of A_infinity
        else
        error('evals wrong')
        end
        
        lambda = next_lambda;
        next_ksi_plus=[ -sqrt(-2*lambda); 0];                                  %stable evec of A_infinity
        
        lambda = last_lambda;
        last_ksi_plus=[ -sqrt(-2*lambda); 0];                                  %stable evec of A_infinity

        %next_lambda... integration of u from x=L -> L_prime
        xspan = [L L_prime];           %x integration range
        lambda = next_lambda;
        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;
        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;
        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)

            V = D*(exp(-2*omega*x) - 2*exp(-omega*x));
            B = 2*(lambda-V);
            A = [	0   -B
                    1    0	];
            u_x = A*u;

        end
        %----------------------------------------------------------
    end

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

end