PRO MODIS_APPLY_LUT, FILE_LIST, BAND_LIST, DATA_TYPE, REF_LIST, LUT

;- Apply destriping LUT to one or more MODIS L1B input files

COMPILE_OPT IDL2

;- Check arguments
if (n_elements(file_list) eq 0) then message, 'Argument FILE_LIST is undefined'
if (n_elements(band_list) eq 0) then message, 'Argument BAND_LIST is undefined'
if (n_elements(data_type) eq 0) then message, 'Argument DATA_TYPE is undefined'
if (n_elements(ref_list) eq 0) then message, 'Argument REF_LIST is undefined'
if (n_elements(lut) eq 0) then message, 'Argument LUT is undefined'

;- Get scale factor for 1KM, HKM, QKM images
case data_type of
  '1KM' : scale = 1
  'HKM' : scale = 2
  'QKM' : scale = 4
endcase

;- Set number of detectors and pixels per scan
ndet = scale * 10
npix = scale * 1354

;- Loop over input files
for file_index = 0, n_elements(file_list) - 1 do begin

  ;- Get file name
  file = file_list[file_index]
  
  ;- Get the mirror side data (dimensions: nrow)
  mirror = read_modis_mirror(file)

  ;- Resize mirror side data (dimensions: nscans)
  nscan = n_elements(mirror)
  nrow = nscan * ndet
  mirror = rebin(mirror, nrow, /sample)
    
  ;- Open the file for HDF SDS read/write access
  hdfid = hdf_sd_start(file, /rdwr)
  
  ;- Create row index array
  row_index = lindgen(nrow)

  ;- Loop over bands
  for band_index = 0, n_elements(band_list) - 1 do begin
    
    ;- Get band number
    band = band_list[band_index]
    
    ;- Get image data (dimensions : npix by nrow)
    read_modis_radiance, hdfid, band, data_type, npix, nrow, image, /raw

    ;- Create output image
    destripe = image
    
    ;- Get reference detector
    ref = ref_list[band_index]

    ;- Loop over mirror sides
    for side = 0, 1 do begin

      ;- Loop over detectors
      for det = 0, 9 do begin

        ;- Get row indices for this mirror side and this detector
        loc = where((mirror eq side) and ((row_index mod ndet) eq det))

        ;- If this detector is not the reference detector, apply destriping LUT
        test = (side eq 0) and (det eq ref)
        if (test eq 0) then begin
          destripe[*, loc] = lut[image[*, loc], det, side, band - 1]
        endif
        
      ;- End of loop over detectors  
      endfor

    ;- End of loop over mirror sides
    endfor

    ;- Replace original invalid values (i.e. missing data)
    loc_invalid = where((image lt 0) or (image gt 32767), count_invalid)
    if (count_invalid gt 0) then destripe[loc_invalid] = image[loc_invalid]

    ;- Set median value of destriped image to median value of original image
    loc_valid = where((image ge 0) and (image le 32767), count_valid)
    median_old = long(median(image[loc_valid]))
    median_new = long(median(destripe[loc_valid]))
    median_del = median_new - median_old
    if (median_del ne 0) then begin
      destripe[loc_valid] = destripe[loc_valid] - median_del
    endif

    ;- Write destriped image data
    write_modis_radiance, hdfid, band, data_type, npix, nrow, destripe, /raw

  ;- End of loop over bands
  endfor
  
  ;- Close the file for HDF SDS access
  hdf_sd_end, hdfid
          
;- End of loop over input files          
endfor

END