;+
; Procedure to display some simple statistics.
;
; Statistics are done on the contents of the promary DC.  One
; can specify the range over which statistics will be computed in
; any of several ways.  To specify a begin and end x-value in the
; currently displayed x-axis units, use the parameters brange and erange.
; If no x-axis range is given, the
; user is prompted to specify the region of interest using the mouse.
;
; @param brange {in}{optional}{type=float} Starting value in x-axis units
; @param erange {in}{optional}{type=float} Ending value in x-axis units
; @keyword full {in}{optional}{type=boolean} Compute stats for full spectrum?
; @keyword ret {out}{optional}{type=structure} Structure containing results
; @examples
; <pre>
;    getrec,1
;    show               ; Use the plotter to set the X-axis to channels
;    stats,0,99         ; Gets stats for first 100 channels
;    show               ; Use the plotter to set the X-axis to GHz
;    stats,1.420,1.421  ; Gets stats
;    stats              ; User clicks plot to determine stats region
;    stats,ret=mystats  ; Return stats in a structure
;    print,mystats.rms
; </pre>
;
; @version $Id: stats.pro,v 1.8 2005/05/30 04:03:43 bgarwood Exp $
;-

pro stats,brange,erange,full=full,ret=ret

    compile_opt idl2
    common gbtplot_common,mystate,xarray

    if !g.line then begin
        if data_valid(!g.s[0]) le 0 then begin
            message,'No data in primary data container',/info
            return
        endif
        if n_elements(full) ne 0 then begin
           bchan = 0
           echan = n_elements(*!g.s[0].data_ptr)-1
        end else if (n_elements(brange) eq 0 or n_elements(erange) eq 0 ) then begin
           print,'Click twice to define stats region'
           a = click()
           gbtoplot,[a.x,a.x],mystate.yrange,color=!white
           brange = a.x
           a = click()
           gbtoplot,[a.x,a.x],mystate.yrange,color=!white
           erange = a.x
           bchan = xtochan(brange)
           echan = xtochan(erange)
        end else begin
           bchan = xtochan(brange)
           echan = xtochan(erange)
        end
        if bchan gt echan then begin
           tmp = bchan & bchan = echan & echan = tmp
        end
        if bchan lt 0 then bchan = 0
        if bchan gt (n_elements(*!g.s[0].data_ptr)-1) then $
           bchan = n_elements(*!g.s[0].data_ptr)-1
        if echan lt 0 then echan = 0
        if echan gt (n_elements(*!g.s[0].data_ptr)-1) then $
           echan = n_elements(*!g.s[0].data_ptr)-1
        if bchan eq echan then begin
           print,'Stats are not supported for a single channel.'
           print,'Data value at index ',bchan,' = ',(*!g.s[0].data_ptr)[bchan]
           return
        end
        data = (*!g.s[0].data_ptr)[bchan:echan]
;     endif else begin
;         if (n_params() eq 2) then begin
; 	        data = (*!g.c[0].data_ptr)[bchan:echan]
;         endif else begin
; 	        data = (*!g.c[0].data_ptr)
;         endelse
    endif else begin
      print,'Statistics not yet supported for Continuum mode.'
      return
    end
    
    mom = moment(data,mdev=mad)
    len = n_elements(data)
    chansize = abs(chantox(bchan)-chantox(bchan+1))
    xmin = chantox(bchan)
    xmax = chantox(echan)
    if xmin gt xmax then begin
      tmp = xmin & xmin = xmax & xmax = tmp
    end
    area = total(data)*chansize
    print,'    Chans    bchan    echan        Xmin        Xmax        Ymin        Ymax'
    print,n_elements(data),bchan,echan,xmin,xmax,min(data),max(data), $
          format='(1x,I8,1x,I8,1x,I8,1x,G11.5,1x,G11.5,1x,G11.5,1x,G11.5)'
    print
    print,'                       Mean      Median         RMS    Variance        Area'   
    print,mom[0],median(data),stddev(data),mom[1],area, $
          format='(16x,G11.5,1x,G11.5,1x,G11.5,1x,G11.5,1x,G11.5)'

     ret = {bchan:bchan,echan:echan,nchan:n_elements(data), $
            xmin:xmin,xmax:xmax,min:min(data), $
            max:max(data),mean:mom[0],median:median(data),rms:stddev(data), $
            variance:mom[1],area:area}
end