;------------------------------------------------------------------------------
;+
; NAME:
;   lambert_lb2pix
;
; PURPOSE:
;   Transform from galactic (l,b) coordinates to (ix,iy) pixel numbers.
;   Latitude runs clockwise from X-axis for NGP, counterclockwise for SGP.
;   This function returns the ZERO-INDEXED pixel position.
;
;   (CRPIX1-CRVAL1, CRPIX2-CRVAL2) define the 1-indexed central pixel
;   location of the pole.
;   For example, if CRPIX1-CRVAL1=512, then the pole is exactly in the middle
;   of 1-indexed pixel number 512; if CRPIX1-CRVAL1=512.5, then the pole falls
;   between pixel numbers 512 and 513.
;
;   If the FRACTIONAL flag is set, then a fractional pixel position is returned.
;
; CALLING SEQUENCE:
;   lambert_lb2pix, lcoord, bcoord, fitshead, xpix, ypix, [ /fractional ]
;
; INPUTS:
;   lcoord:     Galactic longitude in degrees.
;   bcoord:     Galactic latitude in degrees.
;   fitshead:   FITS header data.
;
; OPTIONAL INPUTS:
;   fractional: Set this to return fractional rather than integer pixel numbers.
;
; OUTPUTS:
;   xpix:       X position in integer pixels in range [0,naxis1-1].
;   ypix:       Y position in integer pixels in range [0,naxis2-1].
;
; PROCEDURES CALLED:
;   sxpar(), lambert_lb2xy
;
; REVISION HISTORY:
;   Written by D. Schlegel, 30 May 1996, Durham
;-
;------------------------------------------------------------------------------
pro lambert_lb2pix, lcoord, bcoord, fitshead, xpix, ypix, $
 fractional=fractional

   ; Need five parameters
   if N_params() LT 5 then begin
      print, 'Syntax - lambert_lb2pix, lcoord, bcoord, fitshead, xpix, ypix'
      return
   endif

   ; Get necessary information from FITS header
   nsgp = sxpar(fitshead, 'LAM_NSGP')
   scale = sxpar(fitshead, 'LAM_SCAL')
   crval1 = sxpar(fitshead, 'CRVAL1')
   crval2 = sxpar(fitshead, 'CRVAL2')
   crpix1 = sxpar(fitshead, 'CRPIX1')
   crpix2 = sxpar(fitshead, 'CRPIX2')

   lambert_lb2xy, lcoord, bcoord, nsgp, scale, xcoord, ycoord

   if (keyword_set(fractional)) then begin
      xpix = xcoord + crpix1 - crval1 - 1.0
      ypix = ycoord + crpix2 - crval2 - 1.0
   endif else begin
      xpix = fix(xcoord + crpix1 - crval1 - 0.5) < fix(crpix1 + scale - 1.5)
      ypix = fix(ycoord + crpix2 - crval2 - 0.5) < fix(crpix2 + scale - 1.5)
   endelse

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