apparent_position Subroutine

public subroutine apparent_position(eph, b_target, et, rv_obs_ssb, r_target, status_ok)

Return the position of a target body relative to an observer, corrected for light time and stellar aberration.

see the SPICELIB routine spkapo (with 'lt+s')

Arguments

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

the ephemeris
type(celestial_body), intent(in) :: b_target

target body
real(kind=wp), intent(in) :: et

observer ephemeris time (sec)
real(kind=wp), intent(in), dimension(6) :: rv_obs_ssb

state of the observer (j2000 frame w.r.t. solar system barycenter)
real(kind=wp), intent(out), dimension(3) :: r_target

apparant state of the target (j2000 frame) Corrected for one-way light time and stellar aberration
logical, intent(out) :: status_ok

true if no problems

Calls

proc~~apparent_position~~CallsGraph proc~apparent_position apparent_position get_r get_r proc~apparent_position->get_r proc~axis_angle_rotation axis_angle_rotation proc~apparent_position->proc~axis_angle_rotation proc~cross cross proc~apparent_position->proc~cross proc~unit unit proc~apparent_position->proc~unit proc~axis_angle_rotation->proc~cross proc~axis_angle_rotation->proc~unit

Called by

proc~~apparent_position~~CalledByGraph proc~apparent_position apparent_position proc~get_sun_fraction get_sun_fraction proc~get_sun_fraction->proc~apparent_position

Source Code

    subroutine apparent_position(eph, b_target, et, rv_obs_ssb, r_target, status_ok)

    class(ephemeris_class),intent(inout) :: eph !! the ephemeris
    type(celestial_body),intent(in) :: b_target !! target body
    real(wp),dimension(6),intent(in) :: rv_obs_ssb !! state of the observer
                                                   !! (j2000 frame w.r.t. solar system barycenter)
    real(wp),intent(in) :: et !! observer ephemeris time (sec)
    real(wp),dimension(3),intent(out) :: r_target !! apparant state of the target (j2000 frame)
                                                  !! Corrected for one-way light time and stellar aberration
    logical,intent(out) :: status_ok !! true if no problems

    real(wp),dimension(3) :: r_targ_ssb !! target body r wrt. ssb
    real(wp) :: lt !! one-way light time [sec]

    ! Find the geometric position of the target body with respect to the
    ! solar system barycenter. Subtract the position of the observer
    ! to get the relative position. Use this to compute the one-way
    ! light time.
    call eph%get_r(et,b_target,body_ssb,r_targ_ssb,status_ok)
    if (.not. status_ok) return
    r_targ_ssb = r_targ_ssb - rv_obs_ssb(1:3) ! relative pos of target
    lt = norm2(r_targ_ssb) / c_light ! light time

    ! To correct for light time, find the position of the target body
    ! at the current epoch minus the one-way light time. Note that
    ! the observer remains where he is.
    call eph%get_r(et-lt,b_target,body_ssb,r_targ_ssb,status_ok)
    if (.not. status_ok) return
    r_targ_ssb = r_targ_ssb - rv_obs_ssb(1:3)

    ! At this point, r_targ_ssb contains the geometric or light-time
    ! corrected position of the target relative to the observer

    ! stellar aberration correction
    r_target = stellar_aberration(r_targ_ssb,rv_obs_ssb(4:6))

    contains

        function stellar_aberration ( pobj, vobs ) result(appobj)
            !!  Correct the apparent position of an object for stellar aberration.
            !!  see SPICELIB routine `STELAB`

            real(wp),dimension(3),intent(in) :: pobj
            real(wp),dimension(3),intent(in) :: vobs
            real(wp),dimension(3) :: appobj

            real(wp),dimension(3) :: u, vbyc,h
            real(wp) :: lensqr, sinphi, phi
            real(wp),parameter :: zero_tol = epsilon(1.0_wp) !! tolerance for zero

            u = unit(pobj)
            vbyc = vobs / c_light
            lensqr = dot_product ( vbyc, vbyc )
            if ( lensqr >= 1.0_wp) error stop 'velocity > speed of light'
            h = cross(u, vbyc)
            sinphi  = norm2 ( h )
            if ( abs(sinphi) > zero_tol ) then  ! if (sinphi /= 0)
                ! rotate the position of the object by phi
                ! radians about h to obtain the apparent position.
                phi = asin ( sinphi )
                call axis_angle_rotation ( pobj, h, phi, appobj )
            else
                ! observer is moving along the line of sight to the object,
                ! and no correction is required
                appobj = pobj
            end if

        end function stellar_aberration

    end subroutine apparent_position