Position-angle from spherical coordinates.
Status: vector/matrix support routine.
The result is the bearing (position angle), in radians, of point B with respect to point A. It is in the range -pi to +pi. The sense is such that if B is a small distance "east" of point A, the bearing is approximately +pi/2.
Zero is returned if the two points are coincident.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | intent(in) | :: | al | longitude of point A (e.g. RA) in radians |
||
| real(kind=wp), | intent(in) | :: | ap | latitude of point A (e.g. Dec) in radians |
||
| real(kind=wp), | intent(in) | :: | bl | longitude of point B |
||
| real(kind=wp), | intent(in) | :: | bp | latitude of point B |
||
| real(kind=wp), | intent(out) | :: | theta | position angle of B with respect to A |
subroutine PAS ( al, ap, bl, bp, theta )
implicit none
real(wp),intent(in) :: al !! longitude of point A (e.g. RA) in radians
real(wp),intent(in) :: ap !! latitude of point A (e.g. Dec) in radians
real(wp),intent(in) :: bl !! longitude of point B
real(wp),intent(in) :: bp !! latitude of point B
real(wp),intent(out) :: theta !! position angle of B with respect to A
real(wp) :: dl, x, y
dl = bl - al
y = sin(dl)*cos(bp)
x = sin(bp)*cos(ap) - cos(bp)*sin(ap)*cos(dl)
if ( x/=0.0_wp .or. y/=0.0_wp ) then
theta = atan2(y,x)
else
theta = 0.0_wp
end if
end subroutine PAS