; example script
; use the loaded snapshots to make a time average density probability
; distributions function (pdf)
; at the end the variable mh holds the function and obin1 the
; corresponding density values.

pro time_averaged_density_pdf
@TOOLS/common_blocks.inc

  this_field = 'Density'
  weight_field = 'volume'
  hi.min1 = 1e-3                ; minimum value
  hi.max1 = 1e3                 ; maximum value
  hi.nbins1 = 100               ; number of bins

  widget_control, base2, sensitive = 0
  ntimes = N_ELEMENTS(times)

  tapdf = fltarr(hi.nbins1)     ; field to store result
  totaltime = 0.

  starttime = 0L
;  if time_index lt ntimes-1 then starttime = time_index

; tells jacques to switch to the histogram tab
  J_event, 'HISTOGRAMS'
  hi.field1 = this_field
  hi.weight = weight_field

  for time_index=starttime,ntimes-1,astep DO BEGIN ; loop over all loaded timesteps
     set_list_widget_values                        ; update gui with the new time
     wait, 1.
     select_current_grids       ;
     check_cancel               ; allow user to type  c in command line to stop 
     if stopit then begin
        stopit = 0
        break
     endif
; tells jacques to draw the histogram
     J_event, 'DRAW_BUTTON'

     totaltime += 1 ; assumes snapshots are equidistant in time
     tapdf += mh

     if verbose then print, 'memory use in Mb', memory()/1024./1e3
  endfor
  time_index = ntimes-1         ; keep time_index in loaded range

  tapdf /= totaltime            ; compute the average

  mh = tapdf                    ;  store in this common block variable
                                ; to make it accesible from
                                ; commandline

; show the plot of the final result
  titlestring1=this_field+' ['+get_unit(this_field,/ulabel)+']'
  plot, obin1, (mh) , xlog=shouldbe_logged(this_field), $
        /ylog, psym=10, $
        xtitle=titlestring1, ytitle='dN', $
        xrange=[min(obin1),max(obin1)], $
        xstyle=1, ystyle=1, $
        charsize=2, thick=2, $
        /noclip, pos=[0.15,0.1,0.95,0.95], $
        ticklen=.02, title='time averaged density pdf'

  time_index = N_elements(times)-1
  widget_control, base2, sensitive = 1

return
end

time_averaged_density_pdf
end