altitude_maintenance Subroutine

    subroutine altitude_maintenance(seg,et0,inc0,ran0,tru0,dt_max,n_dvs,dv_total,xf)

    implicit none

    class(segment),intent(inout)      :: seg
    real(wp),intent(in)               :: et0      !! initial ephemeris time (sec)
    real(wp),intent(in)               :: inc0     !! initial inclination - IAU_MOON of date (deg)
    real(wp),intent(in)               :: ran0     !! initial RAAN - IAU_MOON of date (deg)
    real(wp),intent(in)               :: tru0     !! initial true anomaly - IAU_MOON of date (deg)
    real(wp),intent(in)               :: dt_max   !! how long to propagate (days)
    integer,intent(out)               :: n_dvs    !! number of DV maneuvers performed
    real(wp),intent(out)              :: dv_total !! total DV (km/s)
    real(wp),dimension(6),intent(out) :: xf       !! final state - inertial frame (km, km/s)

    real(wp),dimension(3) :: r !! position vector (km)
    real(wp),dimension(3) :: v !! velocity vector (km/s)
    real(wp),dimension(3,3) :: rotmat  !! rotation matrix from ICRF to IAU_MOON
    real(wp),dimension(6) :: x  !! J2000-Moon state vector
    integer :: idid             !! integrator status flag
    real(wp) :: sma             !! circular orbit semi-major axis (km)
    real(wp) :: min_altitude    !! minimum altitude (km)
    real(wp) :: dv              !! periapsis raise maneuver magnitude (km/s)
    real(wp) :: t               !! integration time (sec from et0)
    real(wp) :: tf              !! final integration time (sec from et0)
    real(wp) :: gval            !! event function value
    real(wp) :: tru             !! true anomaly (deg)

    ! initialize
    n_dvs = 0
    dv_total = zero

    ! get initial state in J2000 - Cartesian for integration
    ! note that inc,ran are in moon-centered-of-date-frame
    ! [circular orbit so p=sma]
    sma = seg%r_moon + seg%nominal_altitude  ! initial orbit sma (circular)
    call orbital_elements_to_rv(body_moon%mu,sma,zero,&
                                inc0*deg2rad,ran0*deg2rad,zero,tru0*deg2rad,r,v)

    ! rotate from body-fixed moon of date to j2000:
    rotmat = icrf_to_iau_moon(et0)   ! rotation matrix from inertial to body-fixed Moon frame
    x(1:3) = matmul(transpose(rotmat),r)
    x(4:6) = matmul(transpose(rotmat),v)  ! because using "of date" iau_moon frame

    ! times are relative to initial epoch (sec)
    seg%et_ref = et0
    t  = zero
    tf = dt_max*day2sec

    ! propagate until the altitude is less than min_altitude
    seg%event = 1

    ! main integration loop:
    call seg%first_call()
    do

        call seg%integrate_to_event(t,x,tf,idid=idid,gval=gval)

        if (idid<0) then

            write(error_unit,'(A,*(I5/))') 'idid: ',idid
            error stop 'error in integrator'

        elseif (idid==2 .or. idid==3) then

            ! if we reached the max time, then we are done, so exit
            exit

        elseif (idid == 1000) then  ! a root has been found

            select case (seg%event)
            case(1)

                ! if we hit the min altitude, then integrate to next apoapsis
                ! and raise the periapsis:
                seg%event = 2 ! propagate until apoapsis
                call seg%first_call()  ! have to restart the integration
                                       ! since we just root solved

                ! compute the maneuver here that will be performed
                ! at next apoapsis to raise the periapsis, no matter what.
                dv = periapsis_raise_maneuver(x,sma)

            case (2)

                ! we have propagated to periapsis or apoapsis, if at apoapsis,
                ! perform a maneuver to raise periapsis
                ! if at periapsis, just continue on.

                ! compute current true anomaly:
                tru = true_anomaly(x)

                if (tru<179.0_wp .or. tru>181.0_wp ) then

                    ! we have to keep integrating, since we stopped at periapsis
                    ! note: this is inefficient since we have to
                    ! restart the integration. it would be better to prevent
                    ! to root solver from stopping here.
                    seg%event = 2  ! continue with mode 2
                    call seg%first_call()  ! have to restart integration
                    cycle

                else ! we are at apoapsis

                    ! perform the maneuver we computed when
                    ! the min altitude was triggered
                    ! apply the maneuver along the current apoapsis velocity vector:
                    x(4:6) = x(4:6) + dv * unit(x(4:6))

                    ! keep track of totals:
                    n_dvs = n_dvs + 1
                    dv_total = dv_total + dv

                    write(output_unit,'(*(A,F12.6))') 'maneuver at ', t*sec2hr, &
                                                      ' hr : TRU = ', tru, &
                                                      ' : dv = ', dv

                end if

                seg%event = 1     ! continue with normal mode
                call seg%first_call()  ! have to restart the integration
                                       ! since we root solved and/or changed
                                       ! the state

            case default
                error stop 'invalid event value in altitude_maintenance'
            end select

        else
            write(error_unit,'(A,I5)') 'unknown exit code from integrator: idid=',idid
            error stop 'error in altitude_maintenance'
        end if

    end do

    ! final state:
    xf = x

    end subroutine altitude_maintenance