pro precess, ra, dec, equinox1, equinox2, PRINT = print, FK4 = FK4
;+
; NAME:
;	PRECESS
; PURPOSE:
;	Precess coordinates from EQUINOX1 to EQUINOX2.  For interactive
;	display, one can use the procedure ASTRO which calls PRECESS or use
;	the /PRINT keyword.   The default (RA,DEC) system is FK5 based on 
;	epoch J2000.0 but FK4 based on B1950.0 is available via the FK4 keyword.
;
; CALLING SEQUENCE:
;	PRECESS, ra, dec, [ equinox1, equinox2, /PRINT, /FK4 ]
;
; INPUT - OUTPUT:
;	RA - Input right ascension in DEGREES (scalar or vector)
;	DEC - Input declination in DEGREES (scalar or vector)
;		NOTE: The input RA and DEC are modified by PRECESS to give the 
;		values after precession.
;
; OPTIONAL INPUTS:
;	EQUINOX1 - Original equinox of coordinates, numeric scalar.  If 
;		omitted, then PRECESS will query for EQUINOX1 and EQUINOX2.
;	EQUINOX2 - Equinox of precessed coordinates.
;
; OPTIONAL INPUT KEYWORDS:
;	PRINT - If this keyword is set and non-zero, then the precessed
;		coordinates are displayed at the terminal.
;       FK4   - If this keyword is set, the FK4 (B1950.0) system
;               will be used otherwise FK5 (J2000.0) will be used instead.
;
; RESTRICTIONS:
;	Accuracy of precession decreases for declination values near 90 
;	degrees.  PRECESS should not be used more than 2.5 centuries from
;	2000 on the FK5 system (1950.0 on the FK4 system).
;
; EXAMPLES:
;	(1) The Pole Star has J2000.0 coordinates (2h, 31m, 46.3s, 
;		89d 15' 50.6"); compute its coordinates at J1985.0
;
;	IDL> precess, ten(2,31,46.3)*15, ten(89,15,50.6), 2000, 1985, /PRINT
;
;		====> 2h 16m 22.73s, 89d 11' 47.3"
;
;	(2) Precess the B1950 coordinates of Eps Ind (RA = 21h 59m,33.053s,
;	DEC = (-56d, 59', 33.053") to equinox B1975.
;
;	IDL> ra = ten(21, 59, 33.053)*15
;	IDL> dec = ten(-56, 59, 33.053)
;	IDL> precess, ra, dec ,1950, 1975, /fk4
;
; PROCEDURE:
;	Algorithm from Computational Spherical Astronomy by Taff (1983), 
;	p. 24. (FK4). FK5 constants from "Astronomical Almanac Explanatory
;       Supplement 1992, page 104 Table 3.211.1.
;
; REVISION HISTORY
;	Written, Wayne Landsman, STI Corporation  August 1986
;	Correct negative output RA values   February 1989
;	Added /PRINT keyword      W. Landsman   November, 1991
;	Cleaned up code a bit   W. Landsman   December 1992
;       Provided FK5 (J2000.0)  I. Freedman   January 1994
;-    
  On_error,2                                           ;Return to caller

  npar = N_params()
  deg_to_rad = 0.17453292519943D-1

   if ( npar LT 2 ) then begin 

     print,'Syntax - precess, ra, dec, [ equinox1, equinox2, PRINT =, FK4 =]
     print,'    NOTE: RA and DEC must be supplied in DEGREES'
     return 

  endif else if (npar LT 4) then $
      read,'Enter original and new equinox of coordinates: ',equinox1,equinox2 

  npts = min( [N_elements(ra), N_elements(dec)] )
  if npts EQ 0 then $  
       message,'ERROR - Input RA and DEC must be vectors or scalars'

  ra_rad = ra*deg_to_rad	  ;Convert to double precision if not already
  dec_rad = dec*deg_to_rad 

  a = cos( dec_rad )  

 CASE npts of   	         ;Is RA a vector or scalar?

   1:    x = [a*cos(ra_rad), a*sin(ra_rad), sin(dec_rad)] ;input direction 
;                                                          cosines
   else: begin  	

         x = dblarr(npts,3)
         x(0,0) = a*cos(ra_rad)
         x(0,1) = a*sin(ra_rad)
         x(0,2) = sin(dec_rad)
         x = transpose(x)
         end

   ENDCASE  
;
   sec_to_rad = deg_to_rad/3600.d0

   t = 0.001d0*( equinox2 - equinox1)

IF NOT KEYWORD_SET(fk4)  THEN BEGIN
;   MESSAGE,/INFORM,'RA-DEC system defaults to FK5 (J2000.0)'
   st = 0.001d0*( equinox1 - 2000.d0)
;  Compute 3 rotation angles
   A = sec_to_rad * T * (23062.181D0 + ST*(139.656D0 +0.0139D0*ST) $
    + T*(30.188D0 - 0.344D0*ST+17.998D0*T))

   B = sec_to_rad * T * T * (79.280D0 + 0.410D0*ST + 0.205D0*T) + A

   C = sec_to_rad * T * (20043.109D0 - ST*(85.33D0 + 0.217D0*ST) $
      + T*(-42.665D0 - 0.217D0*ST -41.833D0*T))
ENDIF ELSE BEGIN
;   MESSAGE,/INFORM,'FK4 (B1950.0) RA-DEC system selected.'
   st = 0.001d0*( equinox1 - 1900.d0)
;  Compute 3 rotation angles

   A = sec_to_rad * T * (23042.53D0 + ST*(139.75D0 +0.06D0*ST) $
    + T*(30.23D0 - 0.27D0*ST+18.0D0*T))

   B = sec_to_rad * T * T * (79.27D0 + 0.66D0*ST + 0.32D0*T) + A

   C = sec_to_rad * T * (20046.85D0 - ST*(85.33D0 + 0.37D0*ST) $
      + T*(-42.67D0 - 0.37D0*ST -41.8D0*T))

ENDELSE
  sina = sin(a) &  sinb = sin(b)  & sinc = sin(c)
  cosa = cos(a) &  cosb = cos(b)  & cosc = cos(c)

  r = dblarr(3,3)
  r(0,0) = [ cosa*cosb*cosc-sina*sinb, sina*cosb+cosa*sinb*cosc,  cosa*sinc]
  r(0,1) = [-cosa*sinb-sina*cosb*cosc, cosa*cosb-sina*sinb*cosc, -sina*sinc]
  r(0,2) = [-cosb*sinc, -sinb*sinc, cosc]

  x2 = r#x 	;rotate to get output direction cosines

 if npts EQ 1 then begin	         ;Scalar

	ra_rad = atan(x2(1),x2(0))
	dec_rad = asin(x2(2))

 endif else begin	         ;Vector     

	ra_rad = dblarr(npts) + atan(x2(1,*),x2(0,*))
	dec_rad = dblarr(npts) + asin(x2(2,*))

 endelse

  ra = ra_rad/deg_to_rad
  ra = ra + (ra LT 0.)*360.D            ;RA between 0 and 360 degrees
  dec = dec_rad/deg_to_rad

  if keyword_set( PRINT ) then $
      print, 'Equinox (' + strtrim(equinox2,2) + '): ',adstring(ra,dec,1)

  return
  end