LD Subroutine

  1. astro_module.f90
  2. astro_module
  3. LD

public subroutine LD(bm, p, q, e, em, dlim, p1)

Apply light deflection by a solar-system body, as part of transforming coordinate direction into natural direction.

Status: support routine.

Notes

  1. The algorithm is based on Expr. (70) in Klioner (2003) and Expr. (7.63) in the Explanatory Supplement (Urban & Seidelmann 2013), with some rearrangement to minimize the effects of machine precision.

  2. The mass parameter BM can, as required, be adjusted in order to allow for such effects as quadrupole field.

  3. The barycentric position of the deflecting body should ideally correspond to the time of closest approach of the light ray to the body.

  4. The deflection limiter parameter DLIM is phi^2/2, where phi is the angular separation (in radians) between source and body at which limiting is applied. As phi shrinks below the chosen threshold, the deflection is artificially reduced, reaching zero for phi = 0.

  5. The returned vector P1 is not normalized, but the consequential departure from unit magnitude is always negligible.

  6. To accumulate total light deflection taking into account the contributions from several bodies, call the present routine for each body in succession, in decreasing order of distance from the observer.

  7. For efficiency, validation is omitted. The supplied vectors must be of unit magnitude, and the deflection limiter non-zero and positive.

References

  • Urban, S. & Seidelmann, P. K. (eds), Explanatory Supplement to the Astronomical Almanac, 3rd ed., University Science Books (2013).

  • Klioner, Sergei A., "A practical relativistic model for micro- arcsecond astrometry in space", Astr. J. 125, 1580-1597 (2003).

History

  • IAU SOFA revision: 2013 September 3

Arguments

TypeIntentOptionalAttributesName
real(kind=wp), intent(in) :: bm

mass of the gravitating body (solar masses)
real(kind=wp), intent(in), dimension(3):: p

direction from observer to source (unit vector)
real(kind=wp), intent(in), dimension(3):: q

direction from body to source (unit vector)
real(kind=wp), intent(in), dimension(3):: e

direction from body to observer (unit vector)
real(kind=wp), intent(in) :: em

distance from body to observer (au)
real(kind=wp), intent(in) :: dlim

deflection limiter (Note 4)
real(kind=wp), intent(out), dimension(3):: p1

observer to deflected source (unit vector)

Calls

proc~~ld~~CallsGraph proc~ld LD proc~pdp PDP proc~ld->proc~pdp proc~pxp PXP proc~ld->proc~pxp

Called by

proc~~ld~~CalledByGraph proc~ld LD proc~ldsun LDSUN proc~ldsun->proc~ld proc~ldn LDN proc~ldn->proc~ld proc~atciq ATCIQ proc~atciq->proc~ldsun proc~aticqn ATICQN proc~aticqn->proc~ldn proc~atciqn ATCIQN proc~atciqn->proc~ldn proc~aticq ATICQ proc~aticq->proc~ldsun proc~atciqz ATCIQZ proc~atciqz->proc~ldsun proc~atci13 ATCI13 proc~atci13->proc~atciq proc~atco13 ATCO13 proc~atco13->proc~atciq proc~atic13 ATIC13 proc~atic13->proc~aticq proc~atoc13 ATOC13 proc~atoc13->proc~aticq

Contents

Source Code

LD

Source Code

    subroutine LD ( bm, p, q, e, em, dlim, p1 )

    implicit none

    real(wp),intent(in) :: bm !! mass of the gravitating body (solar masses)
    real(wp),dimension(3),intent(in) :: p !! direction from observer to source (unit vector)
    real(wp),dimension(3),intent(in) :: q !! direction from body to source (unit vector)
    real(wp),dimension(3),intent(in) :: e !! direction from body to observer (unit vector)
    real(wp),intent(in) :: em !! distance from body to observer (au)
    real(wp),intent(in) :: dlim !! deflection limiter (Note 4)
    real(wp),dimension(3),intent(out) :: p1 !! observer to deflected source (unit vector)

    !  Schwarzschild radius of the Sun (au)
    !  = 2 * 1.32712440041 D20 / (2.99792458 D8)^2 / 1.49597870700 D11
    real(wp),parameter :: srs = 1.97412574336e-08_wp

    integer :: i
    real(wp) :: qpe(3), qdqpe, w, eq(3), peq(3)

    !  Q . (Q + E).
    do i=1,3
       qpe(i) = q(i) + e(i)
    end do
    call PDP ( q, qpe, qdqpe )

    !  2 x G x BM / ( EM x c^2 x ( Q . (Q + E) ) ).
    w = bm * srs / em / max ( qdqpe, dlim )

    !  P x (E x Q).
    call PXP ( e, q, eq )
    call PXP ( p, eq, peq )

    !  Apply the deflection.
    do i = 1,3
       p1(i) = p(i) + w*peq(i)
    end do

    end subroutine LD