PAP Subroutine

  1. astro_module.f90
  2. astro_module
  3. PAP

public subroutine PAP(a, b, theta)

Position-angle from two p-vectors.

Status: vector/matrix support routine.

Notes

  1. The result is the position angle, in radians, of direction B with respect to direction A. It is in the range -pi to +pi. The sense is such that if B is a small distance "north" of A the position angle is approximately zero, and if B is a small distance "east" of A the position angle is approximately +pi/2.

  2. A and B need not be unit vectors.

  3. Zero is returned if the two directions are the same or if either vector is null.

  4. If A is at a pole, the result is ill-defined.

History

  • IAU SOFA revision: 2006 November 13

Arguments

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

direction of reference point
real(kind=wp), intent(in), dimension(3):: b

direction of point whose PA is required
real(kind=wp), intent(out) :: theta

position angle of B with respect to A (radians)

Calls

proc~~pap~~CallsGraph proc~pap PAP proc~pm PM proc~pap->proc~pm proc~pdp PDP proc~pap->proc~pdp proc~pmp PMP proc~pap->proc~pmp proc~pxp PXP proc~pap->proc~pxp

Contents

Source Code

PAP

Source Code

    subroutine PAP ( a, b, theta )

    implicit none

    real(wp),dimension(3),intent(in) :: a !! direction of reference point
    real(wp),dimension(3),intent(in) :: b !! direction of point whose PA is required
    real(wp),intent(out) :: theta !! position angle of B with respect to A (radians)

    real(wp) :: am, au(3), bm, st, ct, xa, ya, za, eta(3), &
                xi(3), a2b(3)

    !  Modulus and direction of the A vector.
    call PN ( a, am, au )

    !  Modulus of the B vector.
    call PM ( b, bm )

    !  Deal with the case of a null vector.
    if ( am==0.0_wp .or. bm==0.0_wp ) then
       st = 0.0_wp
       ct = 1.0_wp
    else

       !  The "north" axis tangential from A (arbitrary length).
       xa = a(1)
       ya = a(2)
       za = a(3)
       eta(1) = - xa * za
       eta(2) = - ya * za
       eta(3) = xa*xa + ya*ya

       !  The "east" axis tangential from A (same length).
       call PXP ( eta, au, xi )

       !  The vector from A to B.
       call PMP ( b, a, a2b )

       !  Resolve into components along the north and east axes.
       call PDP ( a2b, xi, st )
       call PDP ( a2b, eta, ct )

       !  Deal with degenerate cases.
       if ( st==0.0_wp .and. ct==0.0_wp ) ct = 1.0_wp

    end if

    !  Position angle.
    theta = atan2(st,ct)

    end subroutine PAP