ballistic_derivs Subroutine

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

Equations of motion for a ballistic segment.

Note

In this example, all the segments have the same equations of motion!

Arguments

Type IntentOptional Attributes Name
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 get_acc get_acc proc~ballistic_derivs->get_acc get_r get_r proc~ballistic_derivs->get_r icrf_to_iau_moon icrf_to_iau_moon proc~ballistic_derivs->icrf_to_iau_moon j2000_to_frame j2000_to_frame proc~ballistic_derivs->j2000_to_frame third_body_gravity third_body_gravity proc~ballistic_derivs->third_body_gravity

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(3) :: r_earth_wrt_moon,r_sun_wrt_moon,r_jupiter_wrt_moon
    real(wp),dimension(3,3) :: rotmat
    real(wp),dimension(3) :: a_geopot  !! central body acc
    real(wp),dimension(3) :: a_earth
    real(wp),dimension(3) :: a_sun
    real(wp),dimension(3) :: a_jupiter
    real(wp),dimension(3) :: a_third_body  !! total third-body acc
    real(wp) :: et !! ephemeris time of `t`
    logical :: status_ok

    select type (me)

    class is (segment)

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

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

        if (me%pointmass_central_body) then
            ! pointmass moon gravity model
            a_geopot = -mu_moon / norm2(r)**3 * r
        else
            ! geopotential gravity:
            ! first get the rotation matrix from inertial to body-fixed Moon frame
            select case (grav_frame)
            case (1); rotmat = icrf_to_iau_moon(et) ! iau_moon
            case (2); rotmat = me%moon_pa%j2000_to_frame(et) ! moon_pa (splined)
            case default; error stop 'invalid grav_frame in ballistic_derivs'
            end select
            rb = matmul(rotmat,r)           ! r in body-fixed frame
            call me%grav%get_acc(rb,grav_n,grav_m,a_geopot)  ! get the acc due to the geopotential
            a_geopot = matmul(transpose(rotmat),a_geopot)    ! convert acc back to inertial frame
        end if

        ! third-body state vectors (wrt the central body, which is the moon in this case):
        ! [inertial frame]
        a_third_body  = zero
        if (me%include_pointmass_earth) then
            call me%eph%get_r(et,body_earth,body_moon,r_earth_wrt_moon,status_ok)
            call me%third_body_gravity(r,r_earth_wrt_moon,mu_earth,a_earth)
            a_third_body = a_third_body + a_earth
        end if
        if (me%include_pointmass_sun) then
            call me%eph%get_r(et,body_sun,body_moon,r_sun_wrt_moon,status_ok)
            call me%third_body_gravity(r,r_sun_wrt_moon,mu_sun,a_sun)
            a_third_body = a_third_body + a_sun
        end if
        if (me%include_pointmass_jupiter) then
            call me%eph%get_r(et,body_jupiter,body_moon,r_jupiter_wrt_moon,status_ok)
            call me%third_body_gravity(r,r_jupiter_wrt_moon,mu_jupiter, a_jupiter)
            a_third_body = a_third_body + a_jupiter
        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