Equations of motion for a ballistic orbit around the moon.
| Type | Intent | Optional | 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 () |
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