generate_patch_points Subroutine

    subroutine generate_patch_points(me, lstar, tstar, period, x0, z0, ydot0, pp)

    implicit none

    class(mission_type),intent(inout) :: me
    real(wp),intent(in) :: lstar
    real(wp),intent(in) :: tstar
    real(wp),intent(in) :: period
    real(wp),intent(in) :: x0
    real(wp),intent(in) :: z0
    real(wp),intent(in) :: ydot0
    type(patch_point),dimension(3), intent(out) :: pp !! periapsis, quarter, apoapsis patch points

    real(wp),dimension(6) :: x_wrt_moon_normalized       ! normalized state wrt to moon
    real(wp),dimension(6) :: x_wrt_moon_unnormalized       ! unnormalized state wrt to moon
    real(wp) :: t_unnormalized ! unnormalized time
    type(ddeabm_class) :: prop  !! integrator
    real(wp) :: mu    !! CRTPB parameter
    real(wp),dimension(:),allocatable :: t_crtbp, x_crtbp,y_crtbp,z_crtbp,vx_crtbp,vy_crtbp,vz_crtbp
    real(wp) :: t    !! normalized time for integration
    real(wp) :: tf   !! final time (normalized) for integration
    real(wp) :: dt   !! time step (normalized) for integration
    real(wp),dimension(6) :: x  !! unnormalized state wrt barycenter, for integration
    integer :: idid  !! ddeabm output flag
    integer :: i     !! counter
    real(wp) :: t_, tf_ !! temp time vars for advancing guess by 1/2 period

    integer,parameter  :: n  = 6 !! number of state variables

    ! preallocate the arrays used in the report function:
    allocate(t_crtbp(0),x_crtbp(0),y_crtbp(0),z_crtbp(0),&
             vx_crtbp(0),vy_crtbp(0),vz_crtbp(0))

    ! first generate the patch points in the normalized system wrt the barycenter.
    ! do this with a numerical integration of the CR3BP equations of motion.

    mu = compute_crtpb_parameter(mu_earth,mu_moon)
    t  = zero          ! start at periapsis
    tf = period / two  ! propagate to apoapsis
    dt = period / four ! time step of 1/4 period
    x = [x0, zero, z0, zero, ydot0, zero] ! initial state: normalized wrt barycenter

    ! integrate and report the points at dt steps (periapsis, 1/4, and apoapsis)
    call prop%initialize(n,maxnum=10000,df=func,&
                         rtol=[me%rtol],atol=[me%atol],&
                         report=report)
    if (.not. me%patch_point_file_is_periapsis) then
        ! advance state by 1/2 period to put x at periapsis
        t_ = t ! make temp copies so the originals aren't changed
        tf_ = tf
        call prop%first_call()
        call prop%integrate(t_,x,tf_,idid=idid,integration_mode=1)
    end if
    call prop%first_call()
    call prop%integrate(t,x,tf,idid=idid,integration_mode=2,tstep=dt)

    ! normalized patch points are now in x_crtbp,y_crtbp,z_crtbp,vx_crtbp,vy_crtbp,vz_crtbp
    ! write(*,*) ''
    ! write(*,*) 'generate_patch_points'
    ! write(*,'(A,*(F10.6,1X))') '  t_crtbp  = ', t_crtbp
    ! write(*,'(A,*(F10.6,1X))') '  x_crtbp  = ', x_crtbp
    ! write(*,'(A,*(F10.6,1X))') '  y_crtbp  = ', y_crtbp
    ! write(*,'(A,*(F10.6,1X))') '  z_crtbp  = ', z_crtbp
    ! write(*,'(A,*(F10.6,1X))') '  vx_crtbp = ', vx_crtbp
    ! write(*,'(A,*(F10.6,1X))') '  vy_crtbp = ', vy_crtbp
    ! write(*,'(A,*(F10.6,1X))') '  vz_crtbp = ', vz_crtbp
    ! write(*,*) ''

    ! transform state to moon-centered and unnormalize
    do i = 1, 3 ! periapsis, quarter, apoapsis

        ! transform state to moon-centered rotating frame

        !            y
        !            ^
        !            |
        !  (M1)------+----------------(M2)  --> x
        !       -mu           1-mu
        !

        x_wrt_moon_normalized = [x_crtbp(i) - (1.0_wp - mu), &
                                 y_crtbp(i), &
                                 z_crtbp(i), &
                                 vx_crtbp(i), &
                                 vy_crtbp(i), &
                                 vz_crtbp(i)]

        ! convert to km, km/s, days (see also unnormalize_variables)
        x_wrt_moon_unnormalized(1:3) = x_wrt_moon_normalized(1:3) * lstar         ! unscale distance
        x_wrt_moon_unnormalized(4:6) = x_wrt_moon_normalized(4:6) * (lstar/tstar) ! unscale velocity
        t_unnormalized = t_crtbp(i) * tstar ! unscale time

        ! results:
        select case (me%N_or_S)
        case ('S')
            ! the data in the file is for the South family
            pp(i) = patch_point(t = t_unnormalized*sec2day,&
                                rv = x_wrt_moon_unnormalized)
        case ('N')
            x_wrt_moon_unnormalized(3) = -x_wrt_moon_unnormalized(3)
            x_wrt_moon_unnormalized(6) = -x_wrt_moon_unnormalized(6)
            pp(i) = patch_point(t = t_unnormalized*sec2day,&
                                rv = x_wrt_moon_unnormalized)
        case default
            error stop 'invalid value for N_or_S: '//trim(me%N_or_S)
        end select

    end do

    contains

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

        call crtbp_derivs(mu,x,xdot)

    end subroutine func

    subroutine report(me,t,x)  !! report function
        implicit none
        class(ddeabm_class),intent(inout)    :: me
        real(wp),intent(in)                  :: t
        real(wp),dimension(:),intent(in)     :: x

        t_crtbp   = [t_crtbp,  t]
        x_crtbp   = [x_crtbp,  x(1)]
        y_crtbp   = [y_crtbp,  x(2)]
        z_crtbp   = [z_crtbp,  x(3)]
        vx_crtbp  = [vx_crtbp, x(4)]
        vy_crtbp  = [vy_crtbp, x(5)]
        vz_crtbp  = [vz_crtbp, x(6)]

    end subroutine report

    end subroutine generate_patch_points