/usr/share/gnudatalanguage/astrolib/jprecess.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 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 | pro jprecess, ra, dec, ra_2000, dec_2000, MU_RADEC = mu_radec, $
PARALLAX = parallax, RAD_VEL = rad_vel, EPOCH = epoch
;+
; NAME:
; JPRECESS
; PURPOSE:
; Precess astronomical coordinates from B1950 to J2000
; EXPLANATION:
; Calculate the mean place of a star at J2000.0 on the FK5 system from the
; mean place at B1950.0 on the FK4 system.
;
; Use BPRECESS for the reverse direction J2000 ==> B1950
; CALLING SEQUENCE:
; jprecess, ra, dec, ra_2000, dec_2000, [ MU_RADEC = , PARALLAX =
; RAD_VEL =, EPOCH = ]
;
; INPUTS:
; RA,DEC - input B1950 right ascension and declination in *degrees*.
; Scalar or vector
;
; OUTPUTS:
; RA_2000, DEC_2000 - the corresponding J2000 right ascension and
; declination in *degrees*. Same number of elements as RA,DEC
; but always double precision.
;
; OPTIONAL INPUT-OUTPUT KEYWORDS
; MU_RADEC - 2xN element double precision vector containing the proper
; motion in seconds of arc per tropical *century* in right
; ascension and declination.
; PARALLAX - N_element vector giving stellar parallax (seconds of arc)
; RAD_VEL - N_element vector giving radial velocity in km/s
;
; The values of MU_RADEC, PARALLAX, and RADVEL will all be modified
; upon output to contain the values of these quantities in the
; J2000 system. Values will also be converted to double precision.
; The parallax and radial velocity will have a very minor influence on
; the J2000 position.
;
; EPOCH - scalar giving epoch of original observations, default 1950.0d
; This keyword value is only used if the MU_RADEC keyword is not set.
; NOTES:
; The algorithm is taken from the Explanatory Supplement to the
; Astronomical Almanac 1992, page 184.
; Also see Aoki et al (1983), A&A, 128,263
;
; JPRECESS distinguishes between the following two cases:
; (1) The proper motion is known and non-zero
; (2) the proper motion is unknown or known to be exactly zero (i.e.
; extragalactic radio sources). In this case, the algorithm
; in Appendix 2 of Aoki et al. (1983) is used to ensure that
; the output proper motion is exactly zero. Better precision
; can be achieved in this case by inputting the EPOCH of the
; original observations.
;
; The error in using the IDL procedure PRECESS for converting between
; B1950 and J2000 can be up to 12", mainly in right ascension. If
; better accuracy than this is needed then JPRECESS should be used.
;
; EXAMPLE:
; The SAO catalogue gives the B1950 position and proper motion for the
; star HD 119288. Find the J2000 position.
;
; RA(1950) = 13h 39m 44.526s Dec(1950) = 8d 38' 28.63''
; Mu(RA) = -.0259 s/yr Mu(Dec) = -.093 ''/yr
;
; IDL> mu_radec = 100D* [ -15D*.0259, -0.093 ]
; IDL> ra = ten(13,39,44.526)*15.D
; IDL> dec = ten(8,38,28.63)
; IDL> jprecess, ra, dec, ra2000, dec2000, mu_radec = mu_radec
; IDL> print, adstring(ra2000, dec2000,2)
; ===> 13h 42m 12.740s +08d 23' 17.69"
;
; RESTRICTIONS:
; "When transferring individual observations, as opposed to catalog mean
; place, the safest method is to tranform the observations back to the
; epoch of the observation, on the FK4 system (or in the system that was
; used to to produce the observed mean place), convert to the FK5 system,
; and transform to the the epoch and equinox of J2000.0" -- from the
; Explanatory Supplement (1992), p. 180
;
; REVISION HISTORY:
; Written, W. Landsman September, 1992
; Corrected a couple of typos in M matrix October, 1992
; Vectorized, W. Landsman February, 1994
; Implement Appendix 2 of Aoki et al. (1983) for case where proper
; motion unknown or exactly zero W. Landsman November, 1994
; Converted to IDL V5.0 W. Landsman September 1997
; Fixed typo in updating proper motion W. Landsman April 1999
; Make sure proper motion is floating point W. Landsman December 2000
; Use V6.0 notation W. Landsman Mar 2011
;-
On_error,2
compile_opt idl2
if N_params() LT 4 then begin
print,'Syntax - JPRECESS, ra,dec, ra_2000, dec_2000, [MU_RADEC ='
print,' PARALLAX = , RAD_VEL = ]'
print,'Input RA and Dec should be given in DEGREES for B1950'
print,'Proper motion, MU_RADEC, (optional) in arc seconds per *century*'
print,'Parallax (optional) in arc seconds'
print,'Radial Velocity (optional) in km/s'
return
endif
N = N_elements( ra )
if N EQ 0 then message,'ERROR - first parameter (RA vector) is undefined'
if ~keyword_set( RAD_VEL) then rad_vel = dblarr(N) else begin
rad_vel = rad_vel*1.
if N_elements( RAD_VEL ) NE N then message, $
'ERROR - RAD_VEL keyword vector must contain ' + strtrim(N,2) + ' values'
endelse
if N_elements( MU_RADEC) GT 0 then begin
if (N_elements( mu_radec) NE 2*N ) then message, $
'ERROR - MU_RADEC keyword (proper motion) be dimensioned (2,' + $
strtrim(N,2) + ')'
mu_radec = mu_radec*1. ;Make sure at least float
endif
if N_elements(epoch) EQ 0 then epoch = 1950.0d0
if N_elements( Parallax) EQ 0 then parallax = dblarr(N) else $
parallax = parallax*1.
radeg = 180.D/!DPI
sec_to_radian = 1./radeg/3600.0d0
M = [ [+0.9999256782D, +0.0111820610D, +0.0048579479D, $
-0.000551D, +0.238514D, -0.435623D ], $
[ -0.0111820611D, +0.9999374784D, -0.0000271474D, $
-0.238565D, -0.002667D, +0.012254D ], $
[ -0.0048579477D, -0.0000271765D, +0.9999881997D , $
+0.435739D, -0.008541D, +0.002117D ], $
[ +0.00000242395018D, +0.00000002710663D, +0.00000001177656D, $
+0.99994704D, +0.01118251D, +0.00485767D ], $
[ -0.00000002710663D, +0.00000242397878D, -0.00000000006582D, $
-0.01118251D, +0.99995883D, -0.00002714D ], $
[ -0.00000001177656D, -0.00000000006587D, 0.00000242410173D, $
-0.00485767D, -0.00002718D, 1.00000956D] ]
A = 1D-6*[ -1.62557D, -0.31919D, -0.13843D] ;in radians
A_dot = 1D-3*[1.244D, -1.579D, -0.660D ] ;in arc seconds per century
if epoch NE 1950.0d then $
A = A + sec_to_radian * A_dot * (epoch - 1950.0D)/100.0d
ra_rad = ra/radeg & dec_rad = dec/radeg
cosra = cos( ra_rad ) & sinra = sin( ra_rad )
cosdec = cos( dec_rad ) & sindec = sin( dec_rad )
ra_2000 = ra*0.
dec_2000 = dec*0.
for i = 0l, N-1 do begin
r0 = [ cosra[i]*cosdec[i], sinra[i]*cosdec[i], sindec[i] ]
if ~keyword_set( MU_RADEC) then begin
mu_a = 0.0d0
mu_d = 0.0d0
endif else begin
if (N_elements( mu_radec) NE 2*N ) then message, $
'ERROR - MU_RADEC keyword (proper motion) must be dimensioned (2,' + $
strtrim(N,2) + ')'
mu_a = mu_radec[ 0, i]
mu_d = mu_radec[ 1, i ]
endelse
r0_dot = [ -mu_a*sinra[i]*cosdec[i] - mu_d*cosra[i]*sindec[i], $ ;Velocity vector
mu_a*cosra[i]*cosdec[i] - mu_d*sinra[i]*sindec[i] , $
mu_d*cosdec[i] ] + 21.095 * rad_vel[i] * parallax[i] * r0
; Remove the effects of the E-terms of aberration to form r1 and r1_dot.
r1 = r0 - A + (total(r0 * A))*r0
r1_dot = r0_dot - A_dot + ( total( r0 * A_dot))*r0
R_1 = [r1, r1_dot]
R = M # R_1
if ~keyword_set(mu_RADEC) then begin
rr = [ R[0], R[1], R[2]]
v = [ R[3],R[4],R[5] ]
t = ((epoch - 1950.0d0) - 50.00021d)/100.0d0
rr1 = rr + sec_to_radian*v*t
x = rr1[0] & y = rr1[1] & Z = rr1[2]
endif else begin
x = R[0] & y = R[1] & Z = R[2]
x_dot = R[3] & y_dot= R[4] & z_dot = R[5]
endelse
r2 = x^2 + y^2 + z^2
rmag = sqrt( r2 )
dec_2000[i] = asin( z / rmag)
ra_2000[i] = atan( y, x)
if keyword_set(mu_RADEC) then begin
mu_radec[0, i] = ( x*y_dot - y*x_dot) / ( x^2 + y^2)
mu_radec[1, i] = ( z_dot* (x^2 + y^2) - z*(x*x_dot + y*y_dot) ) / $
( r2*sqrt( x^2 + y^2) )
endif
if parallax[i] GT 0. then begin
rad_vel[i] = ( x*x_dot + y*y_dot + z*z_dot )/ (21.095*Parallax[i]*rmag)
parallax[i] = parallax[i] / rmag
endif
endfor
neg = where( ra_2000 LT 0, NNeg )
if Nneg GT 0 then ra_2000[neg] = ra_2000[neg] + 2.D*!DPI
ra_2000 = ra_2000*radeg & dec_2000 = dec_2000*radeg
; Make output scalar if input was scalar
sz = size(ra)
if sz[0] EQ 0 then begin
ra_2000 = ra_2000[0] & dec_2000 = dec_2000[0]
endif
return
end
|