ballistic_derivs Subroutine

private subroutine ballistic_derivs(me, t, x, xdot)

Equations of motion for a ballistic orbit around the moon.

Arguments

TypeIntentOptionalAttributesName
class(ddeabm_class), intent(inout) :: me
real(kind=wp), intent(in) :: t

time [sec from epoch]

real(kind=wp), intent(in), dimension(:):: x

state [r,v] in inertial frame (moon-centered)

real(kind=wp), intent(out), dimension(:):: xdot

derivative of state ()


Calls

proc~~ballistic_derivs~~CallsGraph proc~ballistic_derivs ballistic_derivs icrf_to_iau_moon icrf_to_iau_moon proc~ballistic_derivs->icrf_to_iau_moon third_body_gravity third_body_gravity proc~ballistic_derivs->third_body_gravity

Contents

Source Code


Source Code

    subroutine ballistic_derivs(me,t,x,xdot)

    implicit none

    class(ddeabm_class),intent(inout) :: me
    real(wp),intent(in)               :: t    !! time [sec from epoch]
    real(wp),dimension(:),intent(in)  :: x    !! state [r,v] in inertial frame (moon-centered)
    real(wp),dimension(:),intent(out) :: xdot !! derivative of state (\( dx/dt \))

    real(wp),dimension(3) :: r,rb,v
    reaL(wp),dimension(6) :: rv_earth_wrt_moon,rv_sun_wrt_moon
    real(wp),dimension(3,3) :: rotmat
    real(wp),dimension(3) :: a_geopot
    real(wp),dimension(3) :: a_earth
    real(wp),dimension(3) :: a_sun
    real(wp),dimension(3) :: a_third_body
    real(wp) :: et !! ephemeris time of `t`
    logical :: status_ok  !! ephemeris status flag

    select type (me)

    class is (segment)

        ! get state:
        r = x(1:3)
        v = x(4:6)

        ! compute ephemeris time [sec]:
        et = t + me%et_ref

        ! geopotential gravity:
        rotmat = icrf_to_iau_moon(et)   ! rotation matrix from inertial to body-fixed Moon frame
        rb = matmul(rotmat,r)           ! r in body-fixed frame
        call me%grav%get_acc(rb,me%grav_n,me%grav_m,a_geopot)  ! get the acc due to the geopotential
        a_geopot = matmul(transpose(rotmat),a_geopot)    ! convert acc back to inertial frame

        if (me%include_third_bodies) then
            ! third-body state vectors (wrt the central body, which is the moon in this case):
            ! [inertial frame]
            call me%eph%get_rv(et,body_earth,body_moon,rv_earth_wrt_moon,status_ok)
            call me%eph%get_rv(et,body_sun,body_moon,rv_sun_wrt_moon,status_ok)

            ! third-body perturbation (earth & sun):
            a_third_body = 0.0_wp
            call third_body_gravity(r,rv_earth_wrt_moon(1:3),body_earth%mu,a_earth)
            call third_body_gravity(r,rv_sun_wrt_moon(1:3),body_sun%mu,a_sun)
            a_third_body = a_earth + a_sun
        else
            a_third_body = zero
        end if

        !total derivative vector:
        xdot(1:3) = v
        xdot(4:6) = a_geopot + a_third_body

    class default
        error stop 'invalid class in ballistic_derivs'
    end select

    end subroutine ballistic_derivs