Viewing contents of file '../idllib/astron/contrib/freudenreich/autohist.pro'
PRO AUTOHIST,V, zx,zy,xx,yy, NOPLOT=whatever
;
;+
; NAME:
;       AUTOHIST
;
; PURPOSE:
;       Draws a histogram using automatic bin-sizing.
; 
; EXPLANATION:
;       AUTOHIST chooses a number of bins (initially, SQRT(2*N). If this leads 
;       to a histogram in which > 1/5 of the central 50% of the bins are empty,
;       it decreases the number of bins and tries again. The minimum # bins 
;       is 5. The max=199.   Called by HISTOGAUSS and HALFAGAUSS.    Note that 
;       the intrinsic HISTOGRAM function is *not* used.
;
; CALLING SEQUENCE:
;       AUTOHIST, Sample, XLines, Ylines, XCenters, YCenters, [/NOPLOT]
;
; INPUT:
;       Sample = the vector to be histogrammed
;
; OUTPUT:
;       XLINES = the x coordinates of the points that trace the rectangular 
;               histogram bins
;       YLINES = the y coordinates. To draw the histogram plot YLINES vs XLINES.
;       XCENTERS = the x values of the bin centers
;       YCENTERS = the corresponding y values
;
; OPTIONAL INPUT KEYWORD:
;       /NOPLOT  If set, nothing is drawn
;
; SUBROUTINE CALLS:
;       MED, which calculates a median
; REVISION HISTORY:
;       Written,   H. Freudenreich, STX, 1/91
;       1998 March 17 - Changed shading of histogram.  RSH, RSTX
;-

 ON_ERROR,2
 minbin = 5

 N = N_ELEMENTS(v)
 nb = FIX(SQRT(2.*N)) < 199
 nb = nb > minbin

 x1 = MIN(V)  &  x2 = MAX(V)

TRYAGAIN:

 yy = FLTARR(nb)
 dx = (x2-x1)/nb
 xx = FindGEN(nb)*dx + dx/2. + x1

 ind = (V-x1)/dx
 Q = WHERE( ind GT (nb-1),Count )   &  IF Count GT 0 THEN ind(Q) = nb-1
 Q = WHERE( ind LT 0,Count )        &  IF Count GT 0 THEN ind(Q) = 0

 FOR I = LONG(0),N-1 DO yy(ind(I)) = yy(ind(I))+1.

; Count the fraction of empty bins in the middle half of the histogram:
 x14 = (xx(nb-1)-xx(0))/4. + x1
 X34 = xx(nb-1) - (xx(nb-1)-xx(0))/4.
 Q = WHERE( (yy EQ 0.) AND (xx GT x14) AND (xx LT X34), Count )
 IF (Count GT nb/10) AND (nb GT MInbIN) THEN BEGIN  ; 20% EMPTY
   nb = 3*nb/4
   IF nb LT (2*N) THEN GOTO, TRYAGAIN
 ENDIF

; Fill in zx,zy:
 mb=2*nb+2
 zx=FLTARR(mb)  &  zy=FLTARR(mb)
 zx(0)=x1       &  zy(0)=0.
 K=0 
 FOR J=0,nb-1 DO BEGIN
  K=K+1  &  zx(K) = xx(J)-dx/2.  &  zy(K)=yy(J)
  K=K+1  &  zx(K) = xx(J)+dx/2.  &  zy(K)=yy(J)
 ENDFOR
 K=K+1    &  zx(K)=x2             &  zy(K)=0.

 IF KEYWORD_SET(WHATEVER) THEN RETURN

; Plot, then fill, the bins:
 ytop = MAX(yy(1:nb-2))
 IF yy(0) GT ytop THEN yy(0) = ytop
 IF yy(nb-1) GT ytop THEN yy(nb-1) = ytop

 PLOT,xx,yy,XRAN=[x1-dx, x2+dx],YRAN=[0.,1.1*ytop],PSYM=10,LINE=0

 FOR J = 0,nb-1 DO BEGIN
 IF yy(J) GT 0 THEN BEGIN
     A = [xx(j)-dx/2.,xx(j)+dx/2.,xx(j)+dx/2.,xx(j)-dx/2.] 
     B = [0.,0.,yy(j),yy(j)]
     POLYFILL,A,B,orientation=45
 ENDIF
 ENDFOR

 RETURN
 END