get_c_cdot_two_body_rotating Subroutine

private subroutine get_c_cdot_two_body_rotating(me, eph, to_icrf, c, cdot, status_ok)

rotation matrix for ROTATING <-> ICRF

Type Bound

two_body_rotating_frame

Arguments

Type IntentOptional Attributes Name
class(two_body_rotating_frame), intent(inout) :: me
class(ephemeris_class), intent(inout) :: eph

for ephemeris computations (assumed to have already been initialized)
logical, intent(in) :: to_icrf
real(kind=wp), intent(out), dimension(3,3) :: c
real(kind=wp), intent(out), optional, dimension(3,3) :: cdot
logical, intent(out) :: status_ok

Calls

proc~~get_c_cdot_two_body_rotating~~CallsGraph proc~get_c_cdot_two_body_rotating transformation_module::two_body_rotating_frame%get_c_cdot_two_body_rotating get_rv get_rv proc~get_c_cdot_two_body_rotating->get_rv proc~cross vector_module::cross proc~get_c_cdot_two_body_rotating->proc~cross proc~cross_matrix vector_module::cross_matrix proc~get_c_cdot_two_body_rotating->proc~cross_matrix proc~unit vector_module::unit proc~get_c_cdot_two_body_rotating->proc~unit

Called by

proc~~get_c_cdot_two_body_rotating~~CalledByGraph proc~get_c_cdot_two_body_rotating transformation_module::two_body_rotating_frame%get_c_cdot_two_body_rotating proc~get_c_cdot_two_body_rotating_pulsating transformation_module::two_body_rotating_pulsating_frame%get_c_cdot_two_body_rotating_pulsating proc~get_c_cdot_two_body_rotating_pulsating->proc~get_c_cdot_two_body_rotating

Source Code

    subroutine get_c_cdot_two_body_rotating(me,eph,to_icrf,c,cdot,status_ok)

    !! rotation matrix for ROTATING <-> ICRF

    implicit none

    class(two_body_rotating_frame),intent(inout) :: me
    class(ephemeris_class),intent(inout)         :: eph     !! for ephemeris computations (assumed to have already been initialized)
    logical,intent(in)                           :: to_icrf
    real(wp),dimension(3,3),intent(out)          :: c
    real(wp),dimension(3,3),intent(out),optional :: cdot
    logical,intent(out)                          :: status_ok

    real(wp),dimension(3) :: r         !! position of secondary body w.r.t. primary body [inertial frame]
    real(wp),dimension(3) :: v         !! velocity of secondary body w.r.t. primary body [inertial frame]
    real(wp),dimension(3) :: h         !! angular momentum vector
    real(wp),dimension(3) :: w         !! angular velocity of frame
    logical               :: need_cdot !! if we need to compute `cdot`
    real(wp),dimension(3) :: rhat      !! `r` unit vector
    real(wp),dimension(3) :: hhat      !! `h` unit vector

    need_cdot = present(cdot)

    ! get position & velocity of secondary body w.r.t. primary body, in the inertial frame
    call eph%get_rv(me%et,me%secondary_body,me%primary_body,me%rv12,status_ok)

    if (status_ok) then

        r      = me%rv12(1:3)
        v      = me%rv12(4:6)
        h      = cross(r,v)
        rhat   = unit(r)
        hhat   = unit(h)
        c(1,:) = rhat
        c(3,:) = hhat
        c(2,:) = cross(hhat,rhat)

        if (need_cdot) then
            w = h / dot_product(r,r)            ! see: https://en.wikipedia.org/wiki/Angular_velocity
            cdot = -matmul(c,cross_matrix(w))   ! see: http://arxiv.org/pdf/1311.6010.pdf
        end if

        if (to_icrf) then
            c = transpose(c)
            if (need_cdot) cdot = transpose(cdot)
        end if

    else
        write(error_unit,'(A)') 'Error in get_c_cdot_two_body_rotating: '//&
                                'Could not compute rotation matrix.'
        c = zero
        if (need_cdot) cdot = zero
    end if

    end subroutine get_c_cdot_two_body_rotating