Viewing contents of file '../idllib/uit/pro/bprecess.pro'
pro Bprecess, ra, dec, ra_1950, dec_1950, MU_RADEC = mu_radec, $
PARALLAX = parallax, RAD_VEL = rad_vel
;+
; NAME
; BPRECESS
; PURPOSE:
; Calculate the mean place of a star at B1950.0 on the FK4 system from the
; mean place at J2000.0 on the FK5 system.
;
; CALLING SEQUENCE:
; bprecess, ra, dec, ra_1950, dec_1950, [ MU_RADEC = , PARALLAX =
; RAD_VEL = ]
;
; INPUTS:
; RA,DEC - Input J2000 right ascension and declination in *degrees*.
; Scalar or N element vector
;
; OUTPUTS:
; RA_1950, DEC_1950 - The corresponding B1950 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
; B1950 system. The parallax and radial velocity will have a very
; minor influence on the B1950 position.
;
; NOTES:
; The algorithm is taken from the Astronomical Almanac 1990, page B42.
; Also see Aoki et al (1983), A&A, 128,263
;
; The error in using the IDL procedure PRECESS for converting between
; B1950 and J1950 can be up to 1.5", mainly in right ascension. If
; better accuracy than this is needed then BPRECESS should be used.
;
; EXAMPLE:
; The SAO2000 catalogue gives the J2000 position and proper motion for
; the star HD 119288. Find the B1950 position.
;
; RA(2000) = 13h 42m 12.740s Dec(2000) = 8d 23' 17.69''
; Mu(RA) = -.0257 s/yr Mu(Dec) = -.090 ''/yr
;
; IDL> mu_radec = 100D* [ -15D*.0257, -0.090 ]
; IDL> ra = ten(13, 42, 12.740)*15.D
; IDL> dec = ten(8, 23, 17.69)
; IDL> bprecess, ra, dec, ra1950, dec1950, mu_radec = mu_radec
; IDL> print, adstring(ra1950, dec1950,2)
; ===> 13h 39m 44.526s +08d 38' 28.63"
;
; REVISION HISTORY:
; Written, W. Landsman October, 1992
; Vectorized, W. Landsman February, 1994
;-
On_error,2
if N_params() LT 4 then begin
print,'Syntax - bprecess, ra,dec, ra_1950, dec_1950, [MU_RADEC =
print,' PARALLAX = , RAD_VEL = ]
print,' Input RA and Dec should be given in DEGREES for J2000'
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 not 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 not keyword_set( MU_RADEC) then mu_radec = dblarr( 2, N) else 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.
endelse
if not keyword_set( Parallax) then parallax = dblarr(N) else $
parallax = parallax*1.
radeg = 180.D/!DPI
M = [ [+0.9999256795D, -0.0111814828D, -0.0048590040D, $
-0.000551D, -0.238560D, +0.435730D ], $
[ +0.0111814828D, +0.9999374849D, -0.0000271557D, $
+0.238509D, -0.002667D, -0.008541D ], $
[ +0.0048590039D, -0.0000271771D, +0.9999881946D , $
-0.435614D, +0.012254D, +0.002117D ], $
[ -0.00000242389840D, +0.00000002710544D, +0.00000001177742D, $
+0.99990432D, -0.01118145D, -0.00485852D ], $
[ -0.00000002710544D, -0.00000242392702D, +0.00000000006585D, $
+0.01118145D, +0.99991613D, -0.00002716D ], $
[ -0.00000001177742D, +0.00000000006585D,-0.00000242404995D, $
+0.00485852D, -0.00002717D, +0.99996684D] ]
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
ra_rad = ra/radeg & dec_rad = dec/radeg
cosra = cos( ra_rad ) & sinra = sin( ra_rad )
cosdec = cos( dec_rad ) & sindec = sin( dec_rad )
dec_1950 = dec*0.
ra_1950 = ra*0.
for i = 0, N-1 do begin
r0 = [ cosra(i)*cosdec(i), sinra(i)*cosdec(i), sindec(i) ]
mu_a = mu_radec( 0, i )
mu_d = mu_radec( 1, i )
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
R_0 = [ r0, r0_dot ]
R_1 = M # R_0
; Include the effects of the E-terms of aberration to form r and r_dot.
x1 = R_1(0) & y1 = R_1(1) & z1 = R_1(2)
r1 = R_1(0:2) & r1_dot = R_1(3:5)
rmag = sqrt( x1^2 + y1^2 + z1^2 )
s1 = r1/rmag & s1_dot = r1_dot/rmag
s = s1
for j = 0,2 do begin
r = s1 + A - (total(s * A))*s
s = r/rmag
endfor
r_dot = s1_dot + A_dot - ( total( s * A_dot))*s
x = r(0) & y = r(1) & z = r(2)
x_dot = r_dot(0) & y_dot= r_dot(1) & z_dot = r_dot(2)
r2 = x^2 + y^2 + z^2
rmag = sqrt( r2 )
dec_1950(i) = asin( z / rmag)
ra_1950(i) = atan( y, x)
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) )
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_1950 LT 0, NNeg )
if Nneg GT 0 then ra_1950(neg) = ra_1950(neg) + 2.D*!DPI
ra_1950 = ra_1950*radeg & dec_1950 = dec_1950*radeg
; Make output scalar if input was scalar
sz = size(ra)
if sz(0) EQ 0 then begin
ra_1950 = ra_1950(0) & dec_1950 = dec_1950(0)
endif
return
end