;-----------------------------------------------------------------------
;+
; NAME:
;   issa_proj_ignom
; PURPOSE:
;   Project points from a gnomic projection to longitude and latitude.
;
;   Inverse gnomic projection used for the IRAS Sky Survey Atlas.
;   Projection formulae taken from the IRAS Explanatory Supplement,
;   Volume 1, page X-31, with a different scale factor.
;
;   Pass the pixel positions (x,y) and the plate center (ra0,dec0)
;   in radians, and the scale factor in pixels/radian.
;
;   For the ISSA 500x500 fields, the scale factor is
;     scale = 40 pixels/degree = 2291.8312 pixels/radian
;   and the passed positions (x,y) are in the range [-249,250].
;
; CALLING SEQUENCE:
;   issa_proj_ignom, x, y, ra0, dec0, scale, ra, dec
;
; INPUTS:
;   x:          X position(s) on map in pixel position from center
;   y:          Y position(s) on map in pixel position from center
;   ra0:        Central RA of projection [radians]
;   dec0:       Central DEC of projection [radians]
;   scale:      Scale factor in pixels/radian
;
; OUTPUTS:
;   ra:         RA in radians
;   dec:        DEC in radians
;
; PROCEDURES CALLED:
;
; REVISION HISTORY:   
;   First written by D. Schlegel, Feb 1993, Berkeley in Fortran 77.
;   Modified for IDL by D. Finkbeiner.
;-
;-----------------------------------------------------------------------
pro issa_proj_ignom, x, y, ra0, dec0, scale, ra, dec

   ; Need 7 parameters
   if N_params() LT 7 then begin
      print, 'Syntax - issa_proj_ignom, x, y, ra0, dec0, scale, ra, dec'
      return
   endif

   ; Flip the sign on y... (why?)
   x = x / scale
   y = -y / scale
   r = sqrt(x*x + y*y)
   angle = atan(r>0)

   yne0 = where(y ne 0)

   ; Use a four-quadrant arc-tangent.  Note atan returns between -pi/2,pi/2
   beta = x-x
   if (yne0(0) ge 0) then begin
      beta(yne0) = atan(-x(yne0) / y(yne0))
      ylt0 = where(y lt 0)
      if (ylt0(0) ge 0) then $
         beta(ylt0) = beta(ylt0) + !dpi
   endif
 
   yeq0 = where(y eq 0)
   if (yeq0(0) ge 0) then $
      beta(yeq0) = !dpi*((x(yeq0) lt 0.) - 0.5d0)

   sdec0 = sin(dec0)
   cdec0 = cos(dec0)
   sangle = sin(angle)
   cangle = cos(angle)
   sbeta = sin(beta)
   cbeta = cos(beta)
   xx = sdec0*sangle*cbeta + cdec0*cangle
   yy = sangle*sbeta

   xxne0 = where(xx ne 0)
   ; Use a four-quadrant arc-tangent.  Note atan returns between -pi/2,pi/2
   if (xxne0(0) ge 0) then begin
      ra=xx-xx
      ra(xxne0) = atan(yy(xxne0) / xx(xxne0))
   endif
   xxlt0 = where(xx lt 0)
   if (xxlt0(0) ge 0) then $
      ra(xxlt0) = ra(xxlt0) + !dpi
 
   xxeq0 = where(xx eq 0)
   if (xxeq0(0) ge 0) then $
      ra(xxeq0) = !dpi*((yy(xxeq0) lt 0) - 0.5d0)

   ra = ra + ra0
   dec = asin(sdec0*cangle - cdec0*sangle*cbeta)

   ; Put RA in the range [0,2*pi]...
   hi = where(ra ge 2*!dpi) 
   if (hi(0) ge 0) then $
      ra(hi) = ra(hi) - 2*!dpi
   lo=where(ra lt 0*!dpi) 
   if (lo(0) ge 0) then $
      ra(lo) = ra(lo) + 2*!dpi

   return
end
;-----------------------------------------------------------------------