

In [ ]:

    
# from http://www.math.rutgers.edu/~jmireles/MatLabCode/EarthSunHaloOrbit_NewtonMethod.m

# This program computes a halo orbit near the Sun/Earth L1 libration point, from an initial 
# guess given in metric units. (The initial guess was provided by the instructor).

# main page: http://www.math.rutgers.edu/~jmireles/matLabPage.html

for i = 1:K
    out_i=i
    
    %We need some infromation from the 'reference trajectory'.  
        t0=0;
        tf=tau_n;
        numSteps=2000;
        tspan=linspace(0, tf, numSteps);
        
        %numerical integration
        y0=x_n';                                          %inital condition
        options=odeset('RelTol',2.5e-14,'AbsTol',1e-22);    %set tolerences
        [t,x_ref]=ode113('CRTBP',tspan,y0,options,[],G,mu);
        
    %Start computing the differential of the Newton condition vector
    
    %Need the STM for nth guess at nth time (tau_n is n+1 th time)   
    
        Phi=stateTransCRTBP(t0, tf, x_n, mu);
        
        %some of the terms depend on the vector field
        f_x=vectorField_CRTBP(x_ref(numSteps,:),G,mu);
        
        %Construct the needed differential
        a11=Phi(4,3);
        a12=Phi(4,5);
        a21=Phi(6,3);
        a22=Phi(6,5);
        a31=Phi(2,3);
        a32=Phi(2,5);
        
        %the differential of the constraint vector
        DF=[a11, a12, f_x(4);
            a21, a22, f_x(6);
            a31, a32, f_x(2)];
        
        %Invert
        D=inv(DF);
        
    %Compute next (or final/target) initial conditions    
     %compute next step
        
  Term_1 = [x_n(3); x_n(5); tau_n];
  Term_2 = -D*[x_ref(numSteps, 4); x_ref(numSteps, 6); x_ref(numSteps, 2)];
       x_star= Term_1 + Term_2;
        
    %set up x_n for next itteration or output
    %Put it all into x_n (which is really x_(n+1) as it occurs at the end
    %if the loop
        x_n=[r0(1);
            0;
            x_star(1);
            0;
            x_star(2);
            0];
        tau_n=x_star(3);
end    

outX=x_n