; literally 'automate'
pro automate

  compute_stats, FILENAME='EXTREME_proj_xy_000100'
  compute_stats, FILENAME='EXTREME_proj_xy_000200'
  compute_stats, FILENAME='EXTREME_proj_xy_000300'

end

; writes the PDF and CDF to an ASCII file
pro write_txt, filename, curlyI, PDF, CDF

  print, 'opening '+filename+' for write...'
  openw, 1, filename
  printf, 1, format='(3A20)', 'curlyI', 'PDF', 'CDF'
  for i = 0, n_elements(curlyI)-1 do begin
     printf, 1, format='(3F20.6)', curlyI[i], PDF[i], CDF[i]
  endfor
  close, 1
  print, '...done writing '+filename

  return

end

; ============
; main pro
; ============
pro compute_stats, FILENAME=filename, DATSETNAME=datasetname, STOP=stop

  ; handle keywords
  if not keyword_set(filename) then filename = 'EXTREME_proj_xy_000100'
  if not keyword_set(datasetname) then datasetname = 'dens_proj_xy'

  ; read HDF5 data
  readhdf5data, filename, datasetname, data

  ; reform array to make it 2D
  data = reform(data)

  ; get the dimensions of the data
  dims = size(data,/dimensions)

  ; log column density contrast
  curlyI = alog ( data / mean(data,/double) )

  ; now compute the first 4 moments
  moment = get_moments(curlyI)
  print, 'Mean               (curly I) = ', moment[0]
  print, 'Standard deviation (curly I) = ', moment[1]
  print, 'Skewnees           (curly I) = ', moment[2]
  print, 'Kurtosis           (curly I) = ', moment[3]

  ; compute the PDF and CDF of curly I
  get_pdf, curlyI, PDF=pdf, CDF=cdf, LOC=loc
  window, /free
  plot, loc, pdf, charsize=2, xtitle='curly I', ytitle='PDF', psym=10

  ; plot CDF
  window, /free
  plot, loc, cdf, charsize=2, xtitle='curly I', ytitle='CDF', psym=10

  ; write data to file
  pdf_filename = filename + '_PDF.dat'
  write_txt, pdf_filename, loc, pdf, cdf

  if keyword_set(stop) then stop

  return

end