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 get_rv get_rv proc~ballistic_derivs->get_rv 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

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) :: r_earth_wrt_moon,r_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

    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:
            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,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]
        associate (eph => me%eph)
            select type (eph)
            type is (jpl_ephemeris)
                call eph%get_rv(et,body_earth,body_moon,rv_earth_wrt_moon,status_ok)
                call eph%get_rv(et,body_sun,body_moon,rv_sun_wrt_moon,status_ok)
                r_earth_wrt_moon = rv_earth_wrt_moon(1:3)
                r_sun_wrt_moon   = rv_sun_wrt_moon(1:3)
            type is (jpl_ephemeris_splined)
                ! for this, we can just get position vector only
                call eph%get_r(et,body_earth,body_moon,r_earth_wrt_moon,status_ok)
                call eph%get_r(et,body_sun,body_moon,r_sun_wrt_moon,status_ok)
            class default
                error stop 'invalid eph class'
            end select
        end associate
        ! third-body perturbation (earth & sun):
        call third_body_gravity(r,r_earth_wrt_moon,mu_earth,a_earth)
        call third_body_gravity(r,r_sun_wrt_moon,mu_sun,a_sun)
        a_third_body = a_earth + a_sun

        !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