SEPP Subroutine

  1. astro_module.f90
  2. astro_module
  3. SEPP

public subroutine SEPP(a, b, s)

Angular separation between two p-vectors.

Status: vector/matrix support routine.

Notes

  1. If either vector is null, a zero result is returned.

  2. The angular separation is most simply formulated in terms of scalar product. However, this gives poor accuracy for angles near zero and pi. The present algorithm uses both cross product and dot product, to deliver full accuracy whatever the size of the angle.

History

  • IAU SOFA revision: 2006 November 13

Arguments

TypeIntentOptionalAttributesName
real(kind=wp), intent(in), dimension(3):: a

first p-vector (not necessarily unit length)
real(kind=wp), intent(in), dimension(3):: b

second p-vector (not necessarily unit length)
real(kind=wp), intent(out) :: s

angular separation (radians, always positive)

Calls

proc~~sepp~~CallsGraph proc~sepp SEPP proc~pxp PXP proc~sepp->proc~pxp proc~pdp PDP proc~sepp->proc~pdp

Called by

proc~~sepp~~CalledByGraph proc~sepp SEPP proc~seps SEPS proc~seps->proc~sepp proc~pmsafe PMSAFE proc~pmsafe->proc~seps

Contents

Source Code

SEPP

Source Code

    subroutine SEPP ( a, b, s )

    implicit none

    real(wp),dimension(3),intent(in) :: a !! first p-vector (not necessarily unit length)
    real(wp),dimension(3),intent(in) :: b !! second p-vector (not necessarily unit length)
    real(wp),intent(out) :: s !! angular separation (radians, always positive)

    real(wp) :: axb(3), ss, cs

    !  Sine of the angle between the vectors, multiplied by the two moduli.
    call PXP ( a, b, axb )
    call PM ( axb, ss )

    !  Cosine of the angle, multiplied by the two moduli.
    call PDP ( a, b, cs )

    !  The angle.
    if ( ss/=0.0_wp .or. cs/=0.0_wp ) then
       s = atan2(ss,cs)
    else
       s = 0.0_wp
    end if

    end subroutine SEPP