halo_to_rv_diffcorr Subroutine

public subroutine halo_to_rv_diffcorr(libpoint, mu1, mu2, dist, A_z, n, t1, rv, info, period)

Uses

  • proc~~halo_to_rv_diffcorr~~UsesGraph proc~halo_to_rv_diffcorr halo_orbit_module::halo_to_rv_diffcorr module~math_module math_module proc~halo_to_rv_diffcorr->module~math_module module~minpack_module minpack_module proc~halo_to_rv_diffcorr->module~minpack_module module~rk_module rk_module proc~halo_to_rv_diffcorr->module~rk_module module~kind_module kind_module module~math_module->module~kind_module module~numbers_module numbers_module module~math_module->module~numbers_module module~minpack_module->module~kind_module module~minpack_module->module~numbers_module module~rk_module->module~kind_module module~rk_module->module~numbers_module iso_fortran_env iso_fortran_env module~kind_module->iso_fortran_env module~numbers_module->module~kind_module

Compute the state vector for a halo orbit. This uses the approximation, which is retargeted in the real CR3BP system to produce a periodic orbit.

Todo

use a variable-step size integrator

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: libpoint

Libration point number: [1,2,3]
real(kind=wp), intent(in) :: mu1

grav param for primary body [km3/s2]
real(kind=wp), intent(in) :: mu2

grav param for secondary body [km3/s2]
real(kind=wp), intent(in) :: dist

distance between bodies [km]
real(kind=wp), intent(in) :: A_z

halo z amplitude [km]
integer, intent(in) :: n

halo family: 1, 3
real(kind=wp), intent(in) :: t1

tau1 [rad]
real(kind=wp), intent(out), dimension(6) :: rv

cr3bp normalized state vector [wrt barycenter]
integer, intent(out) :: info

status code (1=no errors)
real(kind=wp), intent(out), optional :: period

period of halo (normalized time units)

Calls

proc~~halo_to_rv_diffcorr~~CallsGraph proc~halo_to_rv_diffcorr halo_orbit_module::halo_to_rv_diffcorr proc~compute_crtpb_parameter crtbp_module::compute_crtpb_parameter proc~halo_to_rv_diffcorr->proc~compute_crtpb_parameter proc~halo_to_rv halo_orbit_module::halo_to_rv proc~halo_to_rv_diffcorr->proc~halo_to_rv proc~hybrd1 minpack_module::hybrd1 proc~halo_to_rv_diffcorr->proc~hybrd1 proc~initialize~2 rk_module::rk_class%initialize proc~halo_to_rv_diffcorr->proc~initialize~2 proc~integrate_to_event~2 rk_module::rk_class%integrate_to_event proc~halo_to_rv_diffcorr->proc~integrate_to_event~2 proc~integrate~2 rk_module::rk_class%integrate proc~halo_to_rv_diffcorr->proc~integrate~2 proc~wrap_angle math_module::wrap_angle proc~halo_to_rv_diffcorr->proc~wrap_angle proc~halo_to_rv->proc~compute_crtpb_parameter proc~compute_libration_points crtbp_module::compute_libration_points proc~halo_to_rv->proc~compute_libration_points proc~hybrd minpack_module::hybrd proc~hybrd1->proc~hybrd proc~destroy~2 rk_module::rk_class%destroy proc~initialize~2->proc~destroy~2 proc~set_function brent_module::brent_class%set_function proc~integrate_to_event~2->proc~set_function proc~zeroin brent_module::brent_class%zeroin proc~integrate_to_event~2->proc~zeroin step step proc~integrate_to_event~2->step proc~integrate~2->step proc~cube_root math_module::cube_root proc~compute_libration_points->proc~cube_root proc~dogleg minpack_module::dogleg proc~hybrd->proc~dogleg proc~dpmpar minpack_module::dpmpar proc~hybrd->proc~dpmpar proc~enorm minpack_module::enorm proc~hybrd->proc~enorm proc~fdjac1 minpack_module::fdjac1 proc~hybrd->proc~fdjac1 proc~qform minpack_module::qform proc~hybrd->proc~qform proc~qrfac minpack_module::qrfac proc~hybrd->proc~qrfac proc~r1mpyq minpack_module::r1mpyq proc~hybrd->proc~r1mpyq proc~r1updt minpack_module::r1updt proc~hybrd->proc~r1updt proc~dogleg->proc~dpmpar proc~dogleg->proc~enorm proc~fdjac1->proc~dpmpar proc~qrfac->proc~dpmpar proc~qrfac->proc~enorm proc~r1updt->proc~dpmpar

Source Code

    subroutine halo_to_rv_diffcorr(libpoint,mu1,mu2,dist,A_z,n,t1,rv,info,period)

    use rk_module
    use minpack_module
    use math_module,     only: wrap_angle

    implicit none

    integer,intent(in)  :: libpoint         !! Libration point number: [1,2,3]
    real(wp),intent(in) :: mu1              !! grav param for primary body [km3/s2]
    real(wp),intent(in) :: mu2              !! grav param for secondary body [km3/s2]
    real(wp),intent(in) :: dist             !! distance between bodies [km]
    real(wp),intent(in) :: A_z              !! halo z amplitude [km]
    integer,intent(in)  :: n                !! halo family: 1, 3
    real(wp),intent(in) :: t1               !! tau1 [rad]
    real(wp),dimension(6),intent(out) :: rv !! cr3bp normalized state vector
                                            !! [wrt barycenter]
    integer,intent(out) :: info             !! status code (1=no errors)
    real(wp),intent(out),optional :: period !! period of halo (normalized time units)

    integer,parameter  :: n_state_vars    = 6         !! number of state variables in the equations of motion
    integer,parameter  :: n_opt_vars      = 2         !! number of variables in the targeting problem
    real(wp),parameter :: t0              = 0.0_wp    !! initial time (normalized) (epoch doesn't matter for cr3bp)
    real(wp),parameter :: tol             = 1.0e-8_wp !! tolerance for event finding
    real(wp),parameter :: xtol            = 1.0e-6_wp !! tolerance for [[hybrd]]
    integer,parameter  :: maxfev          = 1000      !! max number of function evaluations for [[hybrd]]
    integer,parameter  :: n_steps_per_rev = 100       !! number of integration steps per orbit rev

    type(rk8_10_class)               :: prop   !! integrator
    real(wp),dimension(n_opt_vars)   :: x_vy0  !! variables in the targeting problem (x0 and vy0)
    real(wp),dimension(n_opt_vars)   :: vx_vzf !! constraints in the targeting problem (vxf and vzf)
    real(wp),dimension(n_state_vars) :: x0     !! halo initial guess from richardson approximation
    real(wp),dimension(n_state_vars) :: xf     !! state after 1/2 rev (to get the period)
    real(wp) :: tf_actual     !! 1/2 period for retargeted orbit (normalized time)
    real(wp) :: actual_period !! actual halo orbit period for retargeted orbit (normalized time)
    real(wp) :: approx_period !! period approximation (normalized time)
    real(wp) :: dt_to_t1      !! time from `t1=0` to input `t1`
    real(wp) :: gf            !! function value after 1/2 rev (y-coordinate)
    real(wp) :: tau           !! `t1` wrapped from \(- \pi \) to \( \pi \)
    real(wp) :: dt            !! time step (normalized)
    real(wp) :: tmax          !! max final time for event finding integration
    real(wp) :: mu            !! CRTBP parameter

    ! compute the CRTBP mu parameter:
    mu = compute_crtpb_parameter(mu1,mu2)

    ! first we get the halo state approximation at tau1=0:
    call halo_to_rv(libpoint,mu1,mu2,dist,A_z,n,zero,x0,approx_period)

    ! for now, fixed number of integration steps per period:
    dt = approx_period / real(n_steps_per_rev,wp)
    tmax = two * approx_period ! should be enough to find the x-z crossing

    ! initialize the integrator:
    call prop%initialize(n_state_vars,func,g=xz_plane_crossing)

    ! now, solve for a halo:
    x_vy0 = [x0(1),x0(5)]  ! x0 and vy0
    call hybrd1(halo_fcn,2,x_vy0,vx_vzf,tol=xtol,info=info)
    if (info==1) then ! solution converged

        ! now have the solution at t1=0:
        x0(1) = x_vy0(1)
        x0(5) = x_vy0(2)

        ! this is the t1 we want:
        tau = wrap_angle(t1)

        ! if we need the period:
        if (present(period) .or. tau/=zero) then
            ! integrate to the first x-axis crossings (one half rev):
            ! [need to check output...]
            call prop%integrate_to_event(t0,x0,dt,tmax,tol,tf_actual,xf,gf)
            actual_period = two * tf_actual ! normalized period
        end if

        ! now we want to propagate to the input tau1
        if (tau==zero) then
            ! already have the solution
            rv = x0
        else
            ! now, integrate from t1=0 to input t1 to get rv:
            dt_to_t1 = actual_period * (tau / twopi)
            call prop%integrate(t0,x0,dt,dt_to_t1,rv)
        end if

        if (present(period)) period = actual_period

    else ! there was an error
        write(error_unit,'(A)') 'Error: the halo targeting problem did not converge.'
        rv = x0
    end if

    !call prop%destroy()

    contains
!*******************************************************************************

    !***************************************************************************
        subroutine halo_fcn(n,xvec,fvec,iflag)
        !! Halo function for [[hybrd1]]

        implicit none

        integer,intent(in)                :: n      !! `n=2` in this case
        real(wp),dimension(n),intent(in)  :: xvec   !! x_vy0
        real(wp),dimension(n),intent(out) :: fvec   !! [vxf,vzf]
        integer,intent(inout)             :: iflag  !! status flag (set negative
                                                    !! to terminate solver)

        real(wp) :: gf
        real(wp),dimension(6) :: x,x1,xf

        x    = x0      ! initial guess state (z is held fixed)
        x(1) = xvec(1) ! x0
        x(5) = xvec(2) ! vy0

        !integrate to the next x-z-plane crossing:
        call prop%integrate_to_event(t0,x,dt,tmax,tol,tf_actual,xf,gf)

        !want x and z-velocity at the x-z-plane crossing to be zero:
        fvec = [xf(4),xf(6)]

        end subroutine halo_fcn
    !***************************************************************************

    !***************************************************************************
        subroutine func(me,t,x,xdot)
        !! CRTBP derivative function
        implicit none
        class(rk_class),intent(inout)        :: me
        real(wp),intent(in)                  :: t
        real(wp),dimension(me%n),intent(in)  :: x
        real(wp),dimension(me%n),intent(out) :: xdot

        call crtbp_derivs(mu,x,xdot)

        end subroutine func
    !**************************************************************************

    !***************************************************************************
        subroutine xz_plane_crossing(me,t,x,g)
        !! x-z-plane crossing event function
        implicit none
        class(rk_class),intent(inout)        :: me
        real(wp),intent(in)                  :: t
        real(wp),dimension(me%n),intent(in)  :: x
        real(wp),intent(out)                 :: g

        g = x(2)  ! y = 0 at x-z-plane crossing

        end subroutine xz_plane_crossing
    !***************************************************************************

    end subroutine halo_to_rv_diffcorr