function lambda_plane_eval_search_HSpulse

clear; clf;
L=10;           %U-integration start
L_prime=-L;
xspan=[L L_prime];  %x integration range
numsteps=100;  %number of steps used in calculating U_lambda
s_numsteps=10000;   %number of steps around soln loop at L_prime
delta = 2*pi/s_numsteps;    %angle increment

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

%setup graphics boxes
subplot(211);
hold on; grid on;
title('\lambda paths and centres in C-plane, r = 0.014'); xlabel Re(\lambda); ylabel Im(\lambda); zlabel \theta(s)2\pi;
subplot(212); title('\theta(s)/2\pi for the various paths'); hold on; grid on; ylabel \theta(s)/2\pi; xlabel s;
f=getframe;

lambda_0_RE_counter=0;

for lambda_0_RE = -1:0.5:0%-0.18:0.02:0.28
    lambda_0_RE_counter = lambda_0_RE_counter+1
    lambda_0_IM_counter = 0;

    for lambda_0_IM = 0.01:0.5:30
        lambda_0_IM_counter = lambda_0_IM_counter+1

        %lambda path definition along loop around lambda_0
        lambda_0 = lambda_0_RE + i*lambda_0_IM; %centre of lambda path circle 
        r=0.3;                                  %radius of lambda path circle

        subplot(211)
        for j=1:numsteps
            lambda_array(j)=lambda_0 + r*exp(2*pi*i*j/numsteps);
        end
        plot3(real(lambda_array),imag(lambda_array),zeros(length(lambda_array)),'r');
        plot3(real(lambda_0), imag(lambda_0),0,'.r');
        

        %lambda path integration
        theta_init=0;           %theta initial value
        sspan=[0 2*pi];
        [s theta] = ode45(@theta_integration, sspan, theta_init);

        subplot(212)
        plot(s,theta/(2*pi))
        delta_theta = (theta(length(theta))-theta(1))/(2*pi);
        subplot(211)
        plot3(real(lambda_0),imag(lambda_0),delta_theta,'*k')
        Z(lambda_0_IM_counter,lambda_0_RE_counter)=delta_theta;
        f=getframe;
    end
end

figure
meshc(Z)


%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