/usr/share/gnudatalanguage/astrolib/xy2ad.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 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 | pro xy2ad, x, y, astr, a, d
;+
; NAME:
; XY2AD
;
; PURPOSE:
; Compute R.A. and Dec from X and Y and a FITS astrometry structure
; EXPLANATION:
; The astrometry structure must first be extracted by EXTAST from a FITS
; header. The offset from the reference pixel is computed and the CD
; matrix is applied. If distortion is present then this is corrected.
; If a WCS projection (Calabretta & Greisen 2002, A&A, 395, 1077) is
; present, then the procedure WCSXY2SPH is used to compute astronomical
; coordinates. Angles are returned in degrees.
;
; XY2AD is meant to be used internal to other procedures.
; For interactive purposes use XYAD.
;
; CALLING SEQUENCE:
; XY2AD, x, y, astr, a, d
;
; INPUTS:
; X - row position in pixels, scalar or vector
; Y - column position in pixels, scalar or vector
; X and Y should be in the standard IDL convention (first pixel is
; 0), and not the FITS convention (first pixel is 1).
; ASTR - astrometry structure, output from EXTAST procedure containing:
; .CD - 2 x 2 array containing the astrometry parameters CD1_1 CD1_2
; in DEGREES/PIXEL CD2_1 CD2_2
; .CDELT - 2 element vector giving physical increment at reference pixel
; .CRPIX - 2 element vector giving X and Y coordinates of reference pixel
; (def = NAXIS/2)
; .CRVAL - 2 element vector giving R.A. and DEC of reference pixel
; in DEGREES
; .CTYPE - 2 element vector giving projection types
; .LONGPOLE - scalar longitude of north pole
; .LATPOLE - scalar giving native latitude of the celestial pole
; .PV2 - Vector of projection parameter associated with latitude axis
; PV2 will have up to 21 elements for the ZPN projection, up to 3
; for the SIN projection and no more than 2 for any other
; projection
;
; Fields added for version 2:
; .PV1 - Vector of projection parameters associated with longitude axis
; .AXES - 2 element integer vector giving the FITS-convention axis
; numbers associated with astrometry, in ascending order.
; Default [1,2].
; .REVERSE - byte, true if first astrometry axis is Dec/latitude
; .COORDSYS - 1 or 2 character code giving coordinate system, including
; 'C' = RA/Dec, 'G' = Galactic, 'E' = Ecliptic, 'X' = unknown.
; .RADECSYS - String giving RA/Dec system e.g. 'FK4', 'ICRS' etc.
; .EQUINOX - Double giving the epoch of the mean equator and equinox
; .DATEOBS - Text string giving (start) date/time of observations
; .MJDOBS - Modified julian date of start of observations.
; .X0Y0 - Implied offset in intermediate world coordinates if user has
; specified a non-standard fiducial point via PV1 and also
; has set PV1_0a =/ 0 to indicate that the offset should be
; applied in order to place CRVAL at the IWC origin.
; Should be *added* to the IWC derived from application of
; CRPIX, CDELT, CD to the pixel coordinates.
;
; .DISTORT - Optional substructure specifying distortion parameters
;
;
; OUTPUT:
; A - R.A. in DEGREES, same number of elements as X and Y
; D - Dec. in DEGREES, same number of elements as X and Y
;
; RESTRICTIONS:
; Note that all angles are in degrees, including CD and CRVAL
; Also note that the CRPIX keyword assumes an FORTRAN type
; array beginning at (1,1), while X and Y give the IDL position
; beginning at (0,0). No parameter checking is performed.
;
; NOTES:
; XY2AD tests for presence of WCS coordinates by the presence of a dash
; in the 5th character position in the value of CTYPE (e.g 'DEC--SIN').
; PROCEDURES USED:
; TAG_EXIST(), WCSXY2SPH, SIP_EVAL(), TPV_EVAL()
; REVISION HISTORY:
; Written by R. Cornett, SASC Tech., 4/7/86
; Converted to IDL by B. Boothman, SASC Tech., 4/21/86
; Perform CD multiplication in degrees W. Landsman Dec 1994
; Understand reversed X,Y (X-Dec, Y-RA) axes, W. Landsman October 1998
; Consistent conversion between CROTA and CD matrix W. Landsman Oct. 2000
; No special case for tangent projection W. Landsman June 2003
; Work for non-WCS coordinate transformations W. Landsman Oct 2004
; Use CRVAL reference point for non-WCS transformation W.L. March 2007
; Use post V6.0 notation W.L. July 2009
; Some optimisation for large input arrays & use of version 2 astr
; structure, J. P. Leahy July 2013
; Evalue TPV distortion (SCAMP) if present W. Landsman Jan 2014
; Support IRAF TNX porjection M. Sullivan U. of Southamptom Mar 2014
; No longer check that CDELT[0] NE 1 W. Landsman Apr 2015
;-
common Broyden_coeff, pv1, pv2 ;Needed for TPV transformation
compile_opt idl2
if N_params() LT 4 then begin
print,'Syntax -- XY2AD, x, y, astr, a, d'
return
endif
Catch, theError
IF theError NE 0 then begin
Catch,/Cancel
void = cgErrorMsg(/quiet)
RETURN
ENDIF
cd = astr.cd
crpix = astr.crpix
cdelt = astr.cdelt
cd[0,0] *= cdelt[0] & cd[0,1] *= cdelt[0]
cd[1,1] *= cdelt[1] & cd[1,0] *= cdelt[1]
xdif = x - (crpix[0]-1)
ydif = y - (crpix[1]-1)
no_PV1 =0 ;Set if PV1 used by TGV distortion
if tag_exist(astr,'DISTORT') && astr.distort.name EQ 'SIP' then begin
distort = astr.distort
a = distort.a
b = distort.b
na = ((size(a,/dimen))[0])
xdif1 = xdif
ydif1 = ydif
for i=0,na-1 do begin
for j=0,na-1 do begin
if a[i,j] NE 0.0 then xdif1 += xdif^i*ydif^j*a[i,j]
if b[i,j] NE 0.0 then ydif1 += xdif^i*ydif^j*b[i,j]
endfor
endfor
xdif = TEMPORARY(xdif1)
ydif = TEMPORARY(ydif1)
ENDIF
astr2 = TAG_EXIST(astr,'AXES') ; version 2 astrometry structure
xsi = cd[0,0]*xdif + cd[0,1]*ydif ;Can't use matrix notation, in
;case X and Y are vectors
eta = cd[1,0]*TEMPORARY(xdif) + cd[1,1]*TEMPORARY(ydif)
if tag_exist(astr,'DISTORT') && astr.distort.name EQ 'TPV' then begin
pv1 = astr.pv1
pv2 = astr.pv2
result = tpv_eval( [[xsi],[eta]])
xsi = reform( result[*,0] )
eta = reform( result[*,1] )
no_PV1 = 1
ctype = strmid(astr.ctype,0,4) + '-TAN'
endif
if tag_exist(astr,'DISTORT') && astr.distort.name EQ 'TNX' then begin
pv1=astr.distort.lngcor
pv2=astr.distort.latcor
result = tnx_eval( [[xsi],[eta]])
xsi = reform( result[*,0] )
eta = reform( result[*,1] )
no_PV1 = 1
ctype = strmid(astr.ctype,0,4) + '-TAN'
endif
if N_elements(ctype) Eq 0 then ctype = astr.ctype
crval = astr.crval
IF astr2 THEN reverse = astr.reverse ELSE BEGIN
coord = strmid(ctype,0,4)
reverse = ((coord[0] EQ 'DEC-') && (coord[1] EQ 'RA--')) || $
((coord[0] EQ 'GLAT') && (coord[1] EQ 'GLON')) || $
((coord[0] EQ 'ELAT') && (coord[1] EQ 'ELON'))
ENDELSE
if reverse then begin
crval = rotate(crval,2)
temp = TEMPORARY(xsi) & xsi = TEMPORARY(eta) & eta = TEMPORARY(temp)
endif
if strmid(ctype[0],4,1) EQ '-' then begin
if no_PV1 then begin ;Set default values for tangent projection
pv1 = [0.0d,0,90.0d,180.0d,90.0d] & pv2 = [0.0d,0.0d]
endif else begin
pv1 = astr.pv1
pv2 = astr.pv2
endelse
if astr2 THEN $
WCSXY2SPH, xsi, eta, a, d, CTYPE = ctype[0:1], PV1 = pv1, $
PV2 = astr.PV2, CRVAL = crval, CRXY = astr.x0y0 $
ELSE $
WCSXY2SPH, xsi, eta, a, d, CTYPE = ctype[0:1], PV2 = pv2, $
LONGPOLE = astr.longpole, CRVAL = crval, LATPOLE = astr.latpole
endif else begin
a = crval[0] + TEMPORARY(xsi) & d = crval[1] + TEMPORARY(eta)
endelse
return
end
|