FUNCTION CONS_DEC,DEC,X,ASTR,ALPHA	  ;Find line of constant Dec
;+
; NAME:
;	CONS_DEC
; PURPOSE:
;	Returns a set of Y pixels values, given an image with astrometry, and
;	either
;	(1)  A set of X pixel values, and a scalar declination value, or
;	(2)  A set of declination values, and a scalar X value
;
;	Form (1) can be used to find the (X,Y) values of a line of constant
;	declination.  Form (2) can be used to find the Y positions of a set
;	declinations, along a line of constant X.
;
; CALLING SEQUENCE:
;	Y = CONS_DEC( DEC, X, CD, ASTR, ALPHA)
;
; INPUTS:
;	DEC - Declination value(s) in DEGREES (-!PI/2 < DEC < !PI/2).  
;		If X is a vector, then DEC must be a scalar.
;	X -   Specified X pixel value(s) for line of constant declination 
;		If DEC is a vector, then X must be a scalar.
;	ASTR - Astrometry structure, as extracted from a FITS header by the
;		procedure EXTAST
; OUTPUT:
;	Y   - Computed set of Y pixel values.  The number of Y values is the
;		same as either DEC or X, whichever is greater.
;
; OPTIONAL OUTPUT:
;	ALPHA - the right ascensions (DEGREES) associated with the (X,Y) points
;
; RESTRICTIONS:
;	Program will have difficulty converging for declination values near
;	90.  For tangent (gnomonic) projections only.  
;
; REVISION HISTORY:
;	Written, Wayne Landsman  STX Co.	                  April 1988
;	Algorithm  adapted from AIPS memo 27. by Eric Griessen
;	Use new astrometry structure,     W. Landsman    HSTX     Jan. 1994
;-
  On_error,2

  if N_params() lt 3 then begin
	print,'Syntax - Y = CONS_DEC( DEC, X, ASTR, [ALPHA] )'
        return, 0
  endif

  RADEG = 180.0D/!DPI

  crpix = astr.crpix
  crval = astr.crval/RADEG
  cdelt = [ [ astr.cdelt(0), 0.],[0., astr.cdelt(1) ] ]
  cdi = invert(cdelt#astr.cd/RADEG)  ;Greisen uses invert of the CD matrix

  xx = x-(crpix(0)-1.)          ;New coordinate origin
  sdel0 = sin(crval(1)) & cdel0 = cos(crval(1))
  a = cdi(0,0)
  b = xx*cdel0 + sdel0*cdi(0,1)
  sign = 2*( a GT 0 ) - 1 

  alpha = crval(0) + atan(b/a) + $ 
     sign*asin( tan(dec/RADEG)* (xx*sdel0-cdel0*cdi(0,1) )/ (sqrt(a^2+b^2)) )
  alpha = alpha * RADEG
  
  if strmid(astr.ctype(0),5,3) EQ 'UIT' then $
  	uit_ad2xy, alpha, dec, astr, x1, y  else $
	    ad2xy, alpha, dec, astr, x1, y

  return,y
  end