;------------------------------------------------------------------------------
;+
; NAME:
;   lambert_getval
;
; PURPOSE:
;   Read one value at a time from NGP+SGP polar projections.
;
;   Bug with FXREAD: if it tries to read the pixel position [0,0],
;   it reads the entire image instead.  This shouldn't matter, since
;   that pixel is usually outside of the valid projection region.
;
; CALLING SEQUENCE:
;   value = lambert_getval(file_n, file_s, lcoord, bcoord, [ path=path, $
;    /interp, /noloop, /verbose ] )
;
; INPUTS:
;   file_n:     File name for NGP
;   file_s:     File name for SGP
;   lcoord:     Galactic longitude(s) in degrees
;   bcoord:     Galactic latitude(s) in degrees
;
; OPTIONAL INPUTS:
;   path:       File name path; default to ''
;   interp:     Set this flag to return a linearly interpolated value
;               from the 4 nearest pixels
;   noloop:     Set this flag to read all values at once without a FOR loop.
;               This is a faster option for reading a large number of values,
;               but requires reading an entire FITS image into memory.
;               (Actually, the smallest possible sub-image is read.)
;   verbose:    Set this flag for verbose output, printing pixel coordinates
;               and map values.  Setting NOLOOP disables this option.
;
; OUTPUTS:
;   value:      Value(s) from Lambert-projection maps.
;
; PROCEDURES CALLED:
;   lambert_lb2pix, djs_fxread
;
; REVISION HISTORY:
;   Written by D. Schlegel, 22 Sep 1997, Durham
;-
;------------------------------------------------------------------------------
function lambert_getval, file_n, file_s, lcoord, bcoord, path=path, $
 interp=interp, noloop=noloop, verbose=verbose

   ; Need 4 parameters
   if N_params() LT 4 then begin
      print, 'Syntax - value = lambert_getval(file_n, file_s, lcoord, bcoord, $'
      print, ' [ path=path, /noloop, /interp ] )'
      return, -1
   endif

   if (NOT keyword_set(path)) then path = ''

   ; Allocate output data array
   if (N_elements(lcoord) EQ 1) then value = 0.0 $
   else value = float(0*lcoord)

   ; Loop through first NGP then SGP
   for iloop=0, 1 do begin

      if (iloop EQ 0) then begin
         indx = where(bcoord GE 0.0)
         fullname = path+file_n
      endif else begin
         indx = where(bcoord LT 0.0)
         fullname = path+file_s
      endelse

      if (indx(0) NE -1) then begin

         ; Read the header only
	 if strmid(fullname,0,6) ne '/mosaic' then $
		fullname = $
		'/mosaic/group/cmbanalysis/forecast/idl/' + $
		fullname
         hdr = headfits(fullname)

         if (NOT keyword_set(interp)) then begin ; NEAREST PIXELS

            ; Determine the nearest pixel coordinates
            lambert_lb2pix, lcoord(indx), bcoord(indx), hdr, xpix, ypix

            if (keyword_set(noloop)) then begin ; READ FULL IMAGE

               xmin = min(xpix)
               xmax = max(xpix)
               ymin = min(ypix)
               ymax = max(ypix)
               djs_fxread, fullname, subimg, hdr, $
                xmin, xmax, ymin, ymax
               value(indx) = subimg(xpix-xmin,ypix-ymin)

            endif else begin ; READ ONE VALUE AT A TIME

               ; Read one pixel value at a time from data file
               for ii=0, N_elements(indx)-1 do begin
                  djs_fxread, fullname, onedat, hdr, $
                   xpix(ii), xpix(ii), ypix(ii), ypix(ii)
                  value(indx(ii)) = onedat
                  if (keyword_set(verbose)) then $
                   print, format='(f8.3,f8.3,i2,i9,i9,e13.5)', $
                    lcoord(indx), bcoord(indx), $
                    iloop, xpix(ii), ypix(ii), value(indx(ii))
               endfor

            endelse

         endif else begin ; INTERPOLATE

            naxis1 = sxpar(hdr, 'NAXIS1')
            naxis2 = sxpar(hdr, 'NAXIS2')

            ; Determine the pixel coordinates for this projection
            lambert_lb2pix, lcoord(indx), bcoord(indx), hdr, xr, yr, $
             /fractional
            xpix1 = fix(xr-0.5)
            ypix1 = fix(yr-0.5)
            dx = xpix1 - xr + 1.5
            dy = ypix1 - yr + 1.5

            ; Force pixel values to fall within the image boundaries.
            ; Any pixels outside the image are changed to the boundary pixels.
            ibad = where(xpix1 LT 0)
            if (ibad(0) NE -1) then begin
               xpix1(ibad) = 0
               dx(ibad) = 1.0
            endif
            ibad = where(ypix1 LT 0)
            if (ibad(0) NE -1) then begin
               ypix1(ibad) = 0
               dy(ibad) = 1.0
            endif
            ibad = where(xpix1 GE naxis1-1)
            if (ibad(0) NE -1) then begin
               xpix1(ibad) = naxis1-2
               dx(ibad) = 0.0
            endif
            ibad = where(ypix1 GE naxis2-1)
            if (ibad(0) NE -1) then begin
               ypix1(ibad) = naxis2-2
               dy(ibad) = 0.0
            endif

            xpix2 = xpix1 + 1
            ypix2 = ypix1 + 1

            ; Create Nx4 arry of weights
            weight = [ [    dx  *    dy  ] , $
                       [ (1-dx) *    dy  ] , $
                       [    dx  * (1-dy) ] , $
                       [ (1-dx) * (1-dy) ] ]

            if (keyword_set(noloop)) then begin ; READ FULL IMAGE

               xmin = min([xpix1,xpix2])
               xmax = max([xpix1,xpix2])
               ymin = min([ypix1,ypix2])
               ymax = max([ypix1,ypix2])
               djs_fxread, fullname, subimg, hdr, $
                xmin, xmax, ymin, ymax
               value(indx) = $
                subimg(xpix1-xmin,ypix1-ymin) * weight(*,0) + $
                subimg(xpix2-xmin,ypix1-ymin) * weight(*,1) + $
                subimg(xpix1-xmin,ypix2-ymin) * weight(*,2) + $
                subimg(xpix2-xmin,ypix2-ymin) * weight(*,3)

            endif else begin ; READ ONE VALUE AT A TIME

               ; Read 2x2 array at a time from data file
               for ii=0, N_elements(indx)-1 do begin
;                  weight = [ [    dx(ii)  *    dy(ii)  ,  $
;                               (1-dx(ii)) *    dy(ii)  ], $
;                             [    dx(ii)  * (1-dy(ii))  , $
;                               (1-dx(ii)) * (1-dy(ii)) ]  ]
                  djs_fxread, fullname, onedat, hdr, $
                   xpix1(ii), xpix2(ii), ypix1(ii), ypix2(ii)
                  value(indx(ii)) = total(onedat(*) * weight(ii,*))

                  if (keyword_set(verbose)) then $
                   print, format='(f8.3,f8.3,i2,f9.3,f9.3,e13.5)', $
                    lcoord(indx), bcoord(indx), $
                    iloop, xr(ii), yr(ii), value(indx(ii))
               endfor

            endelse

         endelse

      endif

   endfor

   return, value
end
;------------------------------------------------------------------------------