function plot_grand_L2, table_file, plot180=plot180, plot360=plot360
;
; Purpose: Read in *.TAB and *.LBL files of a GRaND L2 map product,
;          plot the data values (optional), and returns a structure 
;          containing grids, values, and image arrays in Claudia and
;          Claudia Double Prime (CDP) coordinate systems. 
;          This function was intended to operate using only 
;          internal IDL routines.
;
; Inputs:
;  table_file - A GRaND L2 table file with an extention of ".TAB". If
;               users provide '' (an empty string) as an input data
;               file, a dialog will pop up to select a file.
;  plot180    - (optional keyword) plots data in Claudia Double Prime
;               coordinate system with east longitude of -180 to 180
;               deg.
;  plot360    - (optional keyword) plots data in Claudia Double Prime
;               coordinate system with east longitude of 0 to 360
;               deg. 
;
; Outputs:
;  results - a structure containing: 
;                grid structure in Claudia coordinate system
;                grid structure in CDP coordinate system
;                data name(s)
;                data values
;                data image array(s) of 720x360 in Claudia coordinate system
;                data image array(s) of 720x360 in CDP coordinate system 
;
; Examples:
;  Download the GRD L2 data from the PDS Small Bodies Node (DWNVGRD_2.zip from
;  http://sbn.psi.edu/pds/resource/dwnvgrd.html & unzip).  The resulting 
;  directory was 'mypath\DWNGRD_2\' (note the use of PC path_separator '\').
;  
;  Example 1 - Fe map
;  
;  ; Read the Fe corrected counts and convert the 720x360 Claudia-Double-Prime map to Fe concentration
;  cd, 'mypath\DWNGRD_2\DATA\'
;  r_Fe=plot_grand_L2('.\IRON\GRD_IRON_CORRECTED_COUNTS_MAP\GRD_IRON_CORRECTED_COUNTS_MAP.TAB', /plot180)
;  ; The structure r_Fe contains a 0.5-deg map of counting rates. I want a map of Fe concentrations.
;  ; The documentation explains how to convert counts/s to concentration (under CATALOG, read
;  ; GRAND_VESTA_IRON_CORR_CNTS_MAP_DS.CAT).  
;  map_Fe=189.2d0*r_Fe.IMAGE0_CDP  ; the result is a 720x360 map of Fe concentrations in wt%.
;  ; Display the data using the IDL image function
;  i=image(map_Fe,rgb_table=39, position=[0.2,0.05,0.9,0.95])
;  c=colorbar(target=i, orientation=1,position=[0.1,0.2,0.15,0.8])
;  
;  Example 2 - Neutron absorption map
;  
;  ; Read and calibrate the thermal neutron absorption map
;  r_abs=plot_grand_L2('.\ABSORPTION\GRD_NEUTRON_ABSORPTION_MAP.TAB', /plot180)
;  ; The structure r_dcp contains 0.5-deg map of the neutron absorption parameter DCP,
;  ; which can be converted to macroscopic absorption units following the 
;  ; instruction in CATALOG, GRAND_VESTA_NEUTRON_ABSORP_MAP.CAT):
;  map_abs=216.4*r_abs.IMAGE0_CDP+66.3  ; (units of 10E-4 square-centimeters/gram)
;  ; Display the data using the IDL image function
;  i=image(map_abs,rgb_table=39, position=[0.2,0.05,0.9,0.95])
;  c=colorbar(target=i, orientation=1,position=[0.1,0.2,0.15,0.8])
;  
;  greatData = plot_grand_L2('myGreatL2Data.TAB', /plot360)
;  greatData = plot_grand_L2('.\IRON\GRD_IRON_CORRECTED_COUNTS_MAP\GRD_IRON_CORRECTED_COUNTS_MAP.TAB', /plot180)
;
; Limitations: 
;          - Requires IDL 8.0 or higher to plot (because of IMAGE function)
;          - Table and label files must reside in the same directory.
;          - Fragile - limited error checking
;
; Acknowledgement: The algorithm used in this code is partially based on the one 
;                  used in grid_image.pro from the IDL Toolkit for Mars Odyssey 
;                  GRS and NS Instruments. Therefore, we acknowledge prior work 
;                  by S. Maurice, D. Lawrence, and O. Gasnault.
;
; Version:
;    Version 0.3 Written by N. Yamashita (NY) on Aug. 2, 2013.
;    Version 1.0 T.H. Prettyman revised code documentation, 5-Mar-2014
;
; Review:
;    5-Mar-2014 - THP tested with IDL 8.2 & PDS GRD L2 map products


; Find a *.TAB file to open
if table_file eq '' then table_file = dialog_pickfile(/read, filter='*.TAB')

; Find a *.LBL file to open
label_file = file_dirname(table_file, /mark_directory) + file_basename(table_file, '.TAB') + '.LBL'

; Read the *.LBL file
openr, lun, label_file, /get_lun
line = '' & name = '' & flag_c = 0

while not eof(lun) do begin

 readf, lun, line &  awk = strsplit(line, /extract)
 if n_elements(awk) eq 3 then begin
  if (awk[0] eq 'OBJECT'     and awk[2] eq 'COLUMN') then flag_c = 1
  if (awk[0] eq 'END_OBJECT' and awk[2] eq 'COLUMN') then flag_c = 0
  if (flag_c eq 1 and awk[0] eq 'NAME' and awk[1] eq '=') then name = [name, awk[2]]
 endif

endwhile

name = name[1:-1]
free_lun, lun

; Read the *.TAB file
t = read_ascii(table_file)
n_column = (size(t.(0)))[1]
n_row    = (size(t.(0)))[2]
if n_column ne n_elements(name) then Message, 'Error reading Label and/or Table files'

n_data = n_column - 7
val = fltarr(n_data, n_row)

; Assuming the following column orders.
pixel_index = reform(t.(0)[0,*])
lat_min     = reform(t.(0)[1,*])
lat_max     = reform(t.(0)[2,*])
delta_lat   = reform(t.(0)[3,*])
lon_min     = reform(t.(0)[4,*])
lon_max     = reform(t.(0)[5,*])
delta_lon   = reform(t.(0)[6,*])
for k=0, n_data-1 do val[k,*] = reform(t.(0)[k+7,*])


res = 2
nx = 360*res
ny = 180*res
ind_x = indgen(nx)
ind_y = indgen(ny)
map = fltarr(n_data, nx, ny)    ; holds equal-angle map with lon. from -180 to 180 deg in CDP
map_c  = fltarr(n_data, nx, ny) ; holds equal-angle map with lon. from -180 to 180 deg in Claudia
map_360= fltarr(n_data, nx, ny) ; holds equal-angle map with lon. from    0 to 360 deg in CDP

for i=0, n_elements(pixel_index)-1 do begin

 ; convert lats and lons to indices
 ;
 ind_ymin = (lat_min[i]+90)*res
 ind_ymax = (lat_min[i]+delta_lat[i]+90)*res-1

 ind_xmin = (lon_min[i]+180)*res
 ind_xmax = (lon_min[i]+delta_lon[i]+180)*res-1

 for j=0, n_data-1 do begin

 ; Fill up the map array with the value.
 ; Roll-over at lon=180deg has been taken into account.
 ;
 if ind_xmax ge nx then begin
  map[ j, ind_xmin:nx-1, ind_ymin:ind_ymax ] = val[j, i]
  map[ j, 0:ind_xmax-nx, ind_ymin:ind_ymax ] = val[j, i]
 endif else begin    
  map[ j, ind_xmin:ind_xmax, ind_ymin:ind_ymax ] = val[j, i]
 endelse

 endfor

endfor

if Keyword_set(plot360) then begin
;; Make lons span from 0 to 360 deg.

 for j=0, n_data-1 do begin
  for i=0, ny-1 do map_360[j,*,i] = shift( map[j,*,i], -180*res )
 endfor
 plotMe = map_360

endif else begin
 if Keyword_set(plot180) then plotMe = map
endelse

if Keyword_set(plot360) OR Keyword_set(plot180) then begin

 for j=0, n_data-1 do begin

  ;; IMAGE function requires IDL 8.0 or higher
  ;
  im = image(reform(plotMe[j,*,*]), rgb_table=39, $
         title=name[7+j]+', Claudia Double Prime Coordinate System')
  c = colorbar(target=im)
  if Keyword_set(plot180) then begin
   t1 = text(0.020, 0.15, '-180$\circ$', /normal, font_size=15)
   t2 = text(0.925, 0.15,  '180$\circ$', /normal, font_size=15)
  endif else begin
   t1 = text(0.035, 0.15,    '0$\circ$', /normal, font_size=15)
   t2 = text(0.925, 0.15,  '360$\circ$', /normal, font_size=15)
  endelse

 endfor

endif

;; convert back to Claudia system
;
lon_min_c = lon_min + 210. & index = where(lon_min_c ge 180.) & lon_min_c[index] -= 360.
lon_max_c = lon_max + 210. & index = where(lon_max_c gt 180.) & lon_max_c[index] -= 360.

g_claudia= {lat_min:lat_min, lat_max:lat_max, delta_lat:delta_lat, lon_min:lon_min_c, lon_max:lon_max_c, delta_lon:delta_lon}
g_cdp    = {lat_min:lat_min, lat_max:lat_max, delta_lat:delta_lat, lon_min:lon_min,   lon_max:lon_max,   delta_lon:delta_lon}


for j=0, n_data-1 do begin
 for i=0, 360-1 do map_c[j,*,i] = shift( map[j,*,i], 210*2 )
endfor

name0 = name[7] & val0 = reform(val[0,*]) & image0_claudia = reform(map_c[0,*,*]) & image0_cdp = reform(map[0,*,*])

if n_data eq 1 then $
 results = {g_claudia:g_claudia, g_cdp:g_cdp, $
            name0:name0, val0:val0, image0_claudia:image0_claudia, image0_cdp:image0_cdp}

if n_data eq 2 then begin
 name1 = name[8] & val1 = reform(val[1,*]) & image1_claudia = reform(map_c[1,*,*]) & image1_cdp = reform(map[1,*,*])
 results = {g_claudia:g_claudia, g_cdp:g_cdp, $
            name0:name0, val0:val0, image0_claudia:image0_claudia, image0_cdp:image0_cdp, $
            name1:name1, val1:val1, image1_claudia:image1_claudia, image1_cdp:image1_cdp}
endif

if n_data ge 3 then begin
 name1 = name[8] & val1 = reform(val[1,*]) & image1_claudia = reform(map_c[1,*,*]) & image1_cdp = reform(map[1,*,*])
 name2 = name[9] & val2 = reform(val[2,*]) & image2_claudia = reform(map_c[2,*,*]) & image2_cdp = reform(map[2,*,*])
 results = {g_claudia:g_claudia, g_cdp:g_cdp, $
            name0:name0, val0:val0, image0_claudia:image0_claudia, image0_cdp:image0_cdp, $
            name1:name1, val1:val1, image1_claudia:image1_claudia, image1_cdp:image1_cdp, $
            name2:name2, val2:val2, image2_claudia:image2_claudia, image2_cdp:image2_cdp}
endif

if n_data gt 4 then print, '*** Only the first 3 values out of' +string(n_data, format='(I2)') $
                           + ' were output. ***'

return, results

END