/usr/share/gnudatalanguage/astrolib/posang.pro is in gdl-astrolib 2018.02.16+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 | PRO POSANG,u,ra1,dc1,ra2,dc2,angle
;+
; NAME:
; POSANG
; PURPOSE:
; Computes rigorous position angle of source 2 relative to source 1
;
; EXPLANATION:
; Computes the rigorous position angle of source 2 (with given RA, Dec)
; using source 1 (with given RA, Dec) as the center.
;
; CALLING SEQUENCE:
; POSANG, U, RA1, DC1, RA2, DC2, ANGLE
;
; INPUTS:
; U -- Describes units of inputs and output:
; 0: everything radians
; 1: RAx in decimal hours, DCx in decimal
; degrees, ANGLE in degrees
; RA1 -- Right ascension of point 1
; DC1 -- Declination of point 1
; RA2 -- Right ascension of point 2
; DC2 -- Declination of point 2
;
; OUTPUTS:
; ANGLE-- Angle of the great circle containing [ra2, dc2] from
; the meridian containing [ra1, dc1], in the sense north
; through east rotating about [ra1, dc1]. See U above
; for units.
;
; PROCEDURE:
; The "four-parts formula" from spherical trig (p. 12 of Smart's
; Spherical Astronomy or p. 12 of Green' Spherical Astronomy).
;
; EXAMPLE:
; For the star 56 Per, the Hipparcos catalog gives a position of
; RA = 66.15593384, Dec = 33.94988843 for component A, and
; RA = 66.15646079, Dec = 33.96100069 for component B. What is the
; position angle of B relative to A?
;
; IDL> RA1 = 66.15593384/15.d & DC1 = 33.95988843
; IDL> RA2 = 66.15646079/15.d & DC2 = 33.96100069
; IDL> posang,1,ra1,dc1,ra2,dc2, ang
; will give the answer of ang = 21.4 degrees
; NOTES:
; (1) If RA1,DC1 are scalars, and RA2,DC2 are vectors, then ANGLE is a
; vector giving the position angle between each element of RA2,DC2 and
; RA1,DC1. Similarly, if RA1,DC1 are vectors, and RA2, DC2 are scalars,
; then DIS is a vector giving the position angle of each element of RA1,
; DC1 and RA2, DC2. If both RA1,DC1 and RA2,DC2 are vectors then ANGLE
; is a vector giving the position angle between each element of RA1,DC1
; and the corresponding element of RA2,DC2. If then vectors are not the
; same length, then excess elements of the longer one will be ignored.
;
; (2) Note that POSANG is not commutative -- the position angle between
; A and B is theta, then the position angle between B and A is 180+theta
; PROCEDURE CALLS:
; ISARRAY()
; HISTORY:
; Modified from GCIRC, R. S. Hill, RSTX, 1 Apr. 1998
; Use V6.0 notation W.L. Mar 2011
;
;-
On_error,2 ;Return to caller
compile_opt idl2
npar = N_params()
IF (npar lt 5) THEN BEGIN
print,'Calling sequence: POSANG,U,RA1,DC1,RA2,DC2,ANGLE'
print,' U = 0 ==> Everything in radians'
print, $
' U = 1 ==> RAx decimal hours, DCx decimal degrees, ANGLE degrees'
RETURN
ENDIF
scalar = (~isarray(ra1) ) && (~isarray(ra2) )
IF scalar THEN BEGIN
IF (ra1 eq ra2) && (dc1 eq dc2) THEN BEGIN
angle = 0.0d0
IF npar eq 5 THEN $
print,'Positions are equal: ', ra1, dc1
RETURN
ENDIF
ENDIF
d2r = !DPI/180.0d0
h2r = !DPI/12.0d0
CASE u OF
0: BEGIN
rarad1 = ra1
rarad2 = ra2
dcrad1 = dc1
dcrad2 = dc2
END
1: BEGIN
rarad1 = ra1*h2r
rarad2 = ra2*h2r
dcrad1 = dc1*d2r
dcrad2 = dc2*d2r
END
ELSE: MESSAGE, $
'U must be 0 for radians or 1 for hours, degrees, arcsec'
ENDCASE
radif = rarad2-rarad1
angle = atan(sin(radif),cos(dcrad1)*tan(dcrad2)-sin(dcrad1)*cos(radif))
IF (u ne 0) THEN angle = angle/d2r
IF (npar eq 5) && (scalar) THEN BEGIN
IF (u ne 0) && (abs(angle) ge 0.1) $
THEN fmt = '(F14.8)' $
ELSE fmt = '(E15.8)'
units = (u ne 0) ? ' degrees' : ' radians'
print,'Position angle of target 2 about target 1 is ' $
+ string(angle,format=fmt) + units
ENDIF
RETURN
END
|