;                  "cindi_hdf_read_sample.pro"
;Sample IDL program to read the CINDI daily HDF files. 
;
;Put a list of the HDF files to read in the file "crfile.lis"


id = ''

farray  = strarr(1500)
farray1 = strarr(1500)
fname   = ' '
ifiles = -1           ; create string arrray of file names to be processed
openr,unit1,'crfile.lis',/get_lun    ; this file contains the names of files to be read.
while not eof(unit1) do begin
  READf,UNIT1,fname
  ifiles = ifiles+1
  farray(ifiles) = fname
  position = strpos(fname,'.hdf')
  farray1(ifiles)= strmid(fname,position-16,1)   ; determine if "r" or "s" file
endwhile
close,unit1
free_lun,unit1

for jj=0,ifiles do begin          ;begin major analysis loop, read and plot 1 file/loop.

  id = farray1(jj)
  
  ;open HDF file and read data
  rank     = long(2)
  maxrank  = long(2)
  dimsizes = lonarr(maxrank)
  fname    = farray(jj)
  filelabel = ' '
  print,'Plotting ',fname
  sdi_id  = hdf_sd_start(fname)
  hdf_sd_fileinfo, sdi_id, datasets, attributes
  print, 'datasets = ', datasets, '    attributes = ', attributes
  if (datasets eq 1) then begin
    sdd_id2 = hdf_sd_select(sdi_id,0)
  endif else begin
    sdd_id2 = hdf_sd_select(sdi_id,2)
  endelse
  hdf_sd_getinfo,sdd_id2,dims=dimsizes
  hdf_sd_getdata,sdd_id2,data
  hdf_sd_endaccess,sdd_id2
  hdf_sd_end,sdi_id
  icount=dimsizes(1)-1


  npts   = n_elements(data(0,*))
  icount = npts-1

;check for any missing data
  date0   = reform(data(0,*))
  quality = where (date0 ne -9999.,npts)
  icount  = npts-1

;Load arrays of individual parameters. Note that I did not bother to use all the available parameters.
;Also note that IDL uses array indices that start at 0, while the documentation uses arrays that start at 1. 
  date    =reform(data(0,quality))
  time    =reform(data(1,quality))
  glat    =reform(data(2,quality))
  glon    =reform(data(3,quality))
  alt     =reform(data(4,quality))
  xpos    =reform(data(5,quality))
  ypos    =reform(data(6,quality))
  zpos    =reform(data(7,quality))
  xvel    =reform(data(8,quality))
  yvel    =reform(data(9,quality))
  zvel    =reform(data(10,quality))
  sza     =reform(data(11,quality))
  lst     =reform(data(12,quality))
  apxalt  =reform(data(13,quality))
  apxlat  =reform(data(14,quality))
  apxlon  =reform(data(15,quality))
  ilat    =reform(data(16,quality))
  magincl =reform(data(17,quality))
  mlat    =reform(data(18,quality))
  mlt     =reform(data(19,quality))
  bgnorth =reform(data(20,quality))
  bgeast  =reform(data(21,quality))
  bgdown  =reform(data(22,quality))
  rpaqual  =reform(data(59,quality))  ;RPA quality flag
  dmqual   =reform(data(60,quality))  ;drift meter quality flag
  vx       =reform(data(61,quality))  ;X-component of ion velocity (s/c coordinates)
  vxsig    =reform(data(62,quality))
  vy       =reform(data(63,quality))  ;Y-component of ion velocity (s/c coordinates)
  vysig    =reform(data(64,quality))
  vz       =reform(data(65,quality))  ;Z-component of ion velocity (s/c coordinates)
  vzsig    =reform(data(66,quality))
  rni      =reform(data(67,quality))  ;RPA derived ion density
  rnisig   =reform(data(68,quality))
  ti       =reform(data(69,quality))  ;ion temperature
  tisig    =reform(data(70,quality))
  am1      =reform(data(71,quality))  ; O+
  frac1    =reform(data(72,quality))
  frac1sig =reform(data(73,quality))
  am2      =reform(data(74,quality))  ; H+
  frac2    =reform(data(75,quality))
  frac2sig =reform(data(76,quality))
  am3      =reform(data(77,quality))  ; He+
  frac3    =reform(data(78,quality))
  frac3sig =reform(data(79,quality))
  am4      =reform(data(80,quality))  ; NO+
  frac4    =reform(data(81,quality))
  frac4sig =reform(data(82,quality))
  am5      =reform(data(83,quality))  ; X+
  frac5    =reform(data(84,quality))
  frac5sig =reform(data(85,quality))
  pot      =reform(data(86,quality))
  potsig   =reform(data(87,quality))
  qualdma  =reform(data(88,quality))  ;quality: drift meter density - a
  dena     =reform(data(89,quality))  ;drift meter density - a
  denasig  =reform(data(90,quality))
  qualdmb  =reform(data(91,quality))  ;quality: drift meter density - b
  denb     =reform(data(92,quality))  ;drift meter density - b
  denbsig  =reform(data(93,quality))
  orbit    =reform(data(94,quality))  ;orbit number
  bgx      =reform(data(95,quality))  ;Bgx (sc coordinates)
  bgy      =reform(data(96,quality))  ;Bgy (sc coordinates)
  bgz      =reform(data(97,quality))  ;Bgz (sc coordinates)
  sra      =reform(data(98,quality))  ;sun-ram angle
  zni      =reform(data(99,quality))  ;RPA 0 Volt ion density
  v_zonal  =reform(data(116,quality)) ;v zonal (magnetic coordinates)
  v_paral  =reform(data(117,quality)) ;v parallel (magnetic coordinates)
  v_merid  =reform(data(118,quality)) ;v meidional (magnetic coordinates)
  
  scvel = sqrt(xvel*xvel+yvel*yvel+zvel*zvel)   ;s/c velocity
  frlight = (frac2 + frac3)*rni                 ; Light Ion concentration
  dens    = rni                                 ;total ion density
  fraco   = frac1                               ;O+ fraction
  frach   = frac2                               ;H+ fraction
  frache  = frac3                               ;He+ fraction
  diplat   =reform( (180./!pi)*atan( 0.5*(bgdown/sqrt(bgnorth^2+bgeast^2)) ) )   ;dip latitude
  
;*******************************************
;Insert your plotting or output code here.
;*******************************************

skipfile:
endfor   ;end major analysis loop.(go back and read next file.)


end