;------------------------------------------------------------------------------
;+
; NAME:
;   lambert_xy_to_lb
; PURPOSE:
;   Transform from (x,y) position from origin to galactic (l,b) coordinates
;   in degrees.
;   Latitude runs clockwise from X-axis for NGP, counterclockwise for SGP.
;
; CALLING SEQUENCE:
;   lambert_xy_to_lb, xcoord, ycoord, nsgp, scale, lcoord, bcoord, $
;    [ lbgood=lbgood ]
;
; INPUTS:
;   prefix:     File name prefix, which will have '_n.fits' and '_s.fits'
;               appended.

;   xcoord:     X position in pixels from the center.
;   ycoord:     Y position in pixels from the center.
;   nsgp:       +1 for NGP projection, -1 for SGP.
;   scale:      Radius of b=0 to b=90 degrees in pixels, usually 512.
;
; OUTPUTS:
;   lcoord:     Galactic longitude in degrees.
;   bcoord:     Galactic latitude in degrees.
;
; OPTIONAL OUTPUTS:
;   lbgood:     Return the following:
;                 0 if rho >= sqrt(2)*scale, i.e. an invalid point
;                 1 if sqrt(2)*scale > rho > scale, i.e. in the far hemisphere
;                 2 if rho <= scale, i.e. in the near hemisphere
;               where rho is the distance in pixels from the pole
;
; PROCEDURES CALLED:
;
; REVISION HISTORY:
;   Written by D. Schlegel, 30 May 1996, Durham
;-
;------------------------------------------------------------------------------
pro lambert_xy_to_lb, xcoord, ycoord, nsgp, scale, lcoord, bcoord, $
 lbgood=lbgood

   BAD_GALL = -1.0
   BAD_GALB = -91.0

; Need six parameters
   if N_params() LT 6 then begin
      print, 'Syntax - lambert_xy_to_lb, xcoord, ycoord, nsgp, scale, $'
      print, ' lcoord, bcoord, [ lbgood=lbgood ])'
      return
   endif

   xvec = xcoord / double(scale)
   yvec = ycoord / double(scale)

   rho2 = xvec*xvec + yvec*yvec
   lbgood = (rho2 LT 2.0)
   lbgood = lbgood + (rho2 LE 1.0)
   rho2 = (lbgood NE 0)*rho2 ; Force rho=0 if radius out of valid range

   lcoord = -nsgp * atan(yvec, xvec)
   xvec = 0
   yvec = 0
   bcoord = asin(nsgp * (1. - rho2))
   rho2 = 0

   ; Convert from radians to degrees
   dradeg = 180.d0 / !dpi
   lcoord = lcoord * dradeg
   bcoord = bcoord * dradeg

   ; Put longitude in the range [0,360).
   lcoord = (lcoord > 0.) + (lcoord LT 0.) * (lcoord + 360.)

   ; If the position is invalid, return BAD_GALL and BAD_GALB
   lcoord = (lbgood EQ 0) * BAD_GALL + (lbgood NE 0) * lcoord
   bcoord = (lbgood EQ 0) * BAD_GALB + (lbgood NE 0) * bcoord

   return
end

;------------------------------------------------------------------------------