docs: def jprecess(ra, dec, mu_radec=None, parallax=None, rad_vel=None, epoch=None):
     """
      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
            Converted to python by A. Drlica-Wagner Feb 2014
     """
     if isinstance(ra, ndarray):
        ra = array(ra)
        dec = array(dec)
     else:
        ra = array([ra0])
        dec = array([dec0])
     n = ra.size
      
     if rad_vel is None:   
        rad_vel = zeros(n,dtype=float)
     else:
        if not isinstance(rad_vel, ndarray):
           rad_vel = array([rad_vel],dtype=float)
        if rad_vel.size != n:   
           raise Exception('ERROR - RAD_VEL keyword vector must be of the same length as RA and DEC')
      
     if (mu_radec is not None):   
        if (array(mu_radec).size != 2 * n):   
           raise Exception('ERROR - MU_RADEC keyword (proper motion) be dimensioned (2,' + strtrim(n, 2) + ')')
        mu_radec = mu_radec * 1.
      
     if parallax is None:   
        parallax = zeros(n,dtype=float)
     else:   
        if not isinstance(parallax, ndarray):
           parallax = array([parallax],dtype=float)
      
     if epoch is None:   
        epoch = 1950.0e0
      
     radeg = 180.e0 / pi
     sec_to_radian = 1/radeg/3600.
     #sec_to_radian = lambda x : deg2rad(x/3600.)
      
     m =  array([ 
        array([+0.9999256782e0, +0.0111820610e0, +0.0048579479e0,  \
                  -0.000551e0,     +0.238514e0,     -0.435623e0      ]), 
        array([ -0.0111820611e0, +0.9999374784e0, -0.0000271474e0,  \
                   -0.238565e0,     -0.002667e0,      +0.012254e0    ]), 
        array([ -0.0048579477e0, -0.0000271765e0, +0.9999881997e0 , \
                   +0.435739e0,      -0.008541e0,      +0.002117e0   ]), 
        array([ +0.00000242395018e0, +0.00000002710663e0, +0.00000001177656e0, \
                   +0.99994704e0,    +0.01118251e0,    +0.00485767e0 ]), 
        array([ -0.00000002710663e0, +0.00000242397878e0, -0.00000000006582e0, \
                   -0.01118251e0,     +0.99995883e0,    -0.00002714e0]), 
        array([ -0.00000001177656e0, -0.00000000006587e0, 0.00000242410173e0, \
                   -0.00485767e0,   -0.00002718e0,     1.00000956e0  ]) ])
      
     a = 1e-6 * array([ -1.62557e0, -0.31919e0, -0.13843e0])        #in radians
     a_dot = 1e-3 * array([1.244e0, -1.579e0, -0.660e0 ])           #in arc seconds per century
      
     ra_rad = deg2rad(ra)
     dec_rad = deg2rad(dec)
     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 in range(n):
        r0 = array([ cosra[i]*cosdec[i], sinra[i]*cosdec[i], sindec[i] ])
      
        if (mu_radec is None):   
           mu_a = 0.
           mu_d = 0.
        else:
           mu_a = mu_radec[ i, 0]
           mu_d = mu_radec[ i, 1]
      
        #Velocity vector
        r0_dot = array([-mu_a*sinra[i]*cosdec[i] - mu_d*cosra[i]*sindec[i], \
                           mu_a*cosra[i]*cosdec[i] - mu_d*sinra[i]*sindec[i] , \
                           mu_d*cosdec[i] ])  + 21.095e0 * rad_vel[i] * parallax[i] * r0
      
        r1 = r0 - a + ((a * r0).sum())*r0
        r1_dot = r0_dot - a_dot + (( a * r0).sum())*r0
      
        r_1 = concatenate((r1, r1_dot))
        r = transpose(dot(transpose(m),transpose(r_1)))
      
        if mu_radec is None:   
           rr = r[0:3] 
           v =  r[3:6] 
           t = ((epoch - 1950.0e0) - 50.00021e0)/100.0
           rr1 = rr + sec_to_radian*v*t
           x = rr1[0]  ; y = rr1[1]  ; z = rr1[2]  
        else:
           x = r[0]  ; y = r[1]  ; z = r[2]  
           x_dot = r[3]  ; y_dot= r[4]  ;  z_dot = r[5]
      
        r2 = x**2 + y**2 + z**2
        rmag = sqrt( r2 )
        dec_2000[i] = arcsin(z / rmag)
        ra_2000[i] = arctan2(y, x)
      
        if mu_radec is not None:
           mu_radec[i, 0] = ( x*y_dot - y*x_dot) / ( x**2 + y**2)
           mu_radec[i, 1] = ( z_dot* (x**2 + y**2) - z*(x*x_dot + y*y_dot) ) /  \
                       ( r2*sqrt( x**2 + y**2) )
      
        if parallax[i] > 0.:
           rad_vel[i] = ( x*x_dot + y*y_dot + z*z_dot )/ (21.095*parallax[i]*rmag)
           parallax[i] = parallax[i] / rmag
      
     neg = (ra_2000 < 0)
     if neg.any() > 0:
        ra_2000[neg] = ra_2000[neg] + 2.0 * pi
        
     ra_2000 = ra_2000*radeg ; dec_2000 = dec_2000*radeg
      
     if ra.size == 1:
        ra_2000 = ra_2000[0]     ; dec_2000 = dec_2000[0]
      
     return ra_2000, dec_2000