; plot_hist ; ; see http://www.physics.emory.edu/~weeks/idl ; for more information and software updates ; ; A user-friendly way to plot roughly normal distributed histograms ; easily with optional gaussian fits. ; The relevant data can be returned with the optional parameter 'data' ; and coff. set number of terms in fit (must be 3 or more) using nterms ; ; original version by John C. Crocker ; slightly revised, Feb '99, Eric R. Weeks; /log added 2-23-99 ; pro plot_hist,d,data,coff,nterms=nterms,fit=fit,xrange=xrange,binsize=binsize,\$ noplot=noplot,oplot=oplot,jitter = jitter,center=center,log=log if keyword_set(fit) then center=1; ERW 2-22-99 if not keyword_set(jitter) then jitter = 0 if keyword_set(nterms) then begin if (nterms lt 3 or nterms gt 6) then begin print,'Resetting nterms to 6. nterms must be 3, 4, 5, or 6! See gaussfit.' nterms = 6 endif endif if not keyword_set(nterms) then nterms = 3 ; don't know why, just gotta do it data = float(d) ; get a little info about the distribution n=n_elements(data) logn=alog10(n) min = min(data,max=max) stdev = stdev(data,mean) min = min([min,mean-(6*stdev)]) max = max([max,mean+(4*stdev)]) if keyword_set(xrange) then begin min = xrange(0) max = xrange(1) w=where(data ge min and data le max,ndata) logn = alog10(ndata) stdev=stdev(data(w)) endif ; calculate a nice binsize if the user doesn't give us one if not keyword_set(binsize) then begin binsize=stdev/logn ; reduce the binsize slightly so that its either 1,2 or 5 times a power of ten logbin=alog10(binsize) fp = (logbin-fix(logbin)) ip = fix(logbin) if (fp lt 0) then begin fp = fp+1 & ip = ip-1 endif coeff = 10^(fp) if (coeff lt 2) then fp = 0 if (coeff ge 2) AND (coeff lt 5) then fp = alog10(2) if (coeff ge 5) then fp = alog10(5) binsize = 10^(float(ip + fp)) endif ; adjust the min value down so that 0 is at the center of a partition minbin = long(min/binsize) if minbin lt 0 then minbin = minbin-1 min = minbin * binsize if (keyword_set(center) and (not keyword_set(xrange))) then xrange=[min,max] ; oh yeah, calculate the histogram hist = histogram(data,binsize=binsize,min=min-(binsize*jitter),max=max+binsize*jitter) ; pad out the hist with a zero so the plot doesn't hang in the air! hist=[hist,[0]] np = n_elements(hist) ; make an 'x' vector for the plot and the fit x = (findgen(np)*binsize) + min - (binsize*jitter) if (not keyword_set(center)) then begin w=where(hist gt 0,nw) minw=w(0) > 1 maxw=w(nw-1) < (n_elements(hist)-2) hist=hist(minw-1:maxw+1) x=x(minw-1:maxw+1) endif ; plot the histogram if not keyword_set(noplot) then begin if keyword_set(oplot) then begin oplot,x,hist,psym=10 endif else begin w=where(hist gt 0) if (keyword_set(xrange)) then begin if (keyword_set(log)) then begin plot,x(w),hist(w),xrange=xrange,psym=4,/ylog endif else begin plot,x,hist,xrange=xrange,psym=10 endelse endif else begin if (keyword_set(log)) then begin plot,x(w),hist(w),psym=4,/ylog endif else begin plot,x,hist,psym=10 endelse endelse endelse endif data = [transpose(x),transpose(hist)] ; do the fit, if desired if keyword_set(fit) then begin ft = gaussfit(x,hist,coff,nterms=nterms) data = [transpose(x),transpose(hist),transpose(ft)] if not keyword_set(noplot) then begin oplot,x,ft,linestyle=3 dy = !y.crange(1)-!y.crange(0) dx = !x.crange(1)-!x.crange(0) yy = !y.crange(0) xx = !x.crange(0) if (keyword_set(log)) then begin xyouts, xx + (dx*0.15), 10^(yy + (dy*.90)), "Height = "+\$ strcompress(string(coff(0),format='(g9.4)')) ,alignment=0.5 xyouts, xx + (dx*0.15), 10^(yy + (dy*.85)), "Offset = "+\$ strcompress(string(coff(1),format='(g9.4)')) ,alignment=0.5 xyouts, xx + (dx*0.15), 10^(yy + (dy*.80)), "Sigma = "+\$ strcompress(string(coff(2),format='(g9.4)')) ,alignment=0.5 endif else begin xyouts, xx + (dx*0.15), yy + (dy*.90), "Height = "+\$ strcompress(string(coff(0),format='(g9.4)')) ,alignment=0.5 xyouts, xx + (dx*0.15), yy + (dy*.85), "Offset = "+\$ strcompress(string(coff(1),format='(g9.4)')) ,alignment=0.5 xyouts, xx + (dx*0.15), yy + (dy*.80), "Sigma = "+\$ strcompress(string(coff(2),format='(g9.4)')) ,alignment=0.5 endelse if (nterms gt 3) then begin if (keyword_set(log)) then begin xyouts, xx + (dx*0.15), 10^(yy + (dy*.75)), "Polynomial:" ,alignment=0.5 xyouts, xx + (dx*0.15), 10^(yy + (dy*.70)), \$ strcompress(string(coff(3),format='(g10.4)')) ,alignment=0.5 endif else begin xyouts, xx + (dx*0.15), yy + (dy*.75), "Polynomial:" ,alignment=0.5 xyouts, xx + (dx*0.15), yy + (dy*.70), \$ strcompress(string(coff(3),format='(g10.4)')) ,alignment=0.5 endelse if (nterms gt 4) then begin if (keyword_set(log)) then begin xyouts, xx + (dx*0.15), 10^(yy + (dy*.65)), "+ "+\$ strcompress(string(coff(4),format='(g9.4)'))+" x" ,alignment=0.5 endif else begin xyouts, xx + (dx*0.15), yy + (dy*.65), "+ "+\$ strcompress(string(coff(4),format='(g9.4)'))+" x" ,alignment=0.5 endelse if (nterms gt 5) then begin if (keyword_set(log)) then begin xyouts, xx + (dx*0.15), 10^(yy + (dy*.60)), "+ "+\$ strcompress(string(coff(5),format='(g10.4)'))+" x^2" ,alignment=0.5 endif else begin xyouts, xx + (dx*0.15), yy + (dy*.60), "+ "+\$ strcompress(string(coff(5),format='(g10.4)'))+" x^2" ,alignment=0.5 endelse endif ; five endif ; four endif ; three endif endif end