function histarray,data,pminval=pminval,pmaxval=pmaxval,binsz=binsz,x=x

;Get rid of points we don't want.

data=data(where(data ne 99999.0))
data=data(where(data eq data))
data=data(where(finite(data) eq 1))
ndata=n_elements(data)
if (ndata eq 0) then begin
	print,'No Valid data to histogram'
	x=0
	return,0
endif
if (ndata gt 5) then begin
	; this loop rejects any data that are more than 20 std.
	; deviations away from the mean of the central 80%
	; of the data.  It takes care of some of the strange results
	; that crop up from screwed up fits.
	data=data(sort(data))
	center=data(long(0.1*ndata):long(0.90*ndata))
	moments=moment(center)
	mn=moments(0)
	disp=sqrt(moments(1))
	data=data(where((data gt mn-20*disp) and (data lt mn+20*disp)))
endif
minval=min(data)
maxval=max(data)
if n_elements(pminval) eq 0 then pminval=minval
if n_elements(pmaxval) eq 0 then pmaxval=maxval
if n_elements(binsz) eq 1 then begin
	hist=histogram(data,binsize=binsz,min=pminval,max=pmaxval)
	bins=n_elements(hist)
	x=findgen(bins)*binsz+pminval
	x=[x(0)-binsz,x,x(bins-1)+binsz]+binsz/2
	hist=[0,hist,0]
	return,hist
endif
if n_elements(binsz) eq 0 then begin
	ndata=n_elements(data)
	maxbins=n_elements(uniq(data(sort(data))))
	binsz=(maxval-minval)/maxbins
endif
tryagain:
hist=histogram(data,binsize=binsz)
bins=n_elements(hist)
x=findgen(bins)*binsz+minval
if (bins le 3) then begin
	x=[x(0)-binsz,x,x(bins-1)+binsz]+binsz/2
	hist=[0,hist,0]
	return,hist
endif
reg=hist(1:bins-2)
up=hist(2:bins-1)
down=hist(0:bins-3)
extrema=((up-reg)/abs(up-reg)+(reg-down)/abs(reg-down))/2
nextrema=n_elements(where(extrema eq 0))
if (nextrema) gt 0.4*float(bins) then begin
 	binsz=binsz*1.2
	goto,tryagain
endif
order=alog10(binsz)
mant=abs(order-ceil(order))
if (mant lt 0.30103001) then binsz=10.^(ceil(order)) else $
binsz=10.^(floor(order))*5
minorder=ceil(alog10(binsz))

minval=floor(minval*10.^(-1*minorder))*(10.^minorder)
hist=histogram(data,binsize=binsz,min=pminval,max=pmaxval)

bins=n_elements(hist)
x=findgen(bins)*binsz+pminval
x=[x(0)-binsz,x,x(bins-1)+binsz]+binsz/2
;while x(0) gt pminval do begin
;	x=[x(0)-binsz,x]
;	hist=[0,hist]
;endwhile
;bins=n_elements(hist)
;while x(bins-1) le pmaxval do begin
;	x=[x,x(bins-1)+binsz]
;	hist=[hist,0]
;	bins=n_elements(hist)
;endwhile
hist=[0,hist,0]
return,hist
end