FUNCTION percentile,array,pctile,flag=flag

; This function returns the data value that lies at the specified
; percentile of the total range of the data in the array.  The
; desired percentile (expressed as a fraction, not a percent) is
; entered in the pctile variable.  If the keyword 'flag' is not
; specified, percentiles less than 0.5 will be approached from
; above, and those greater than 0.5 from below.  These directions
; of approach can be overridden by specifying flag=-1 to approach
; from above or flag=1 to approach from below.

; do some error checking
dim=SIZE(array)
IF dim(0) EQ 0 THEN BEGIN
    MESSAGE,'array must be at least 1-d',/INFORMATIONAL
    RETURN,0
    ENDIF
IF (MIN(pctile) LE 0) OR (MAX(pctile) GE 1) THEN BEGIN
    MESSAGE,'All elements of pctile must be in the range 0 < pctile < 1',/INFORMATIONAL
    RETURN,0
    ENDIF
IF KEYWORD_SET(flag) THEN $
    IF (ABS(flag) NE 1) THEN BEGIN
        MESSAGE,'If set, flag must be -1 or 1',/INFORMATIONAL
        RETURN,0
        ENDIF
IF (MIN(array,MAX=mx) EQ mx) THEN BEGIN
    MESSAGE,'array contains only one value: '+STRTRIM(mx,2),/INFORMATIONAL
    IF N_ELEMENTS(pctile) GT 1 $
        THEN RETURN,REPLICATE(mx,N_ELEMENTS(pctile)) $
	ELSE RETURN,mx
    ENDIF

; choose an appropriate binning factor for the histogram
hbin=2^CEIL(ALOG(MAX(array)/32767.)/ALOG(2.))>1

h=HISTOGRAM(array,MIN=0,BIN=hbin)	; histogram of supplied array
ch=h  

; next line computes cumulative pixel value pdf
FOR jj=1,N_ELEMENTS(ch)-1 DO ch(jj)=ch(jj)+ch(jj-1)

nch=ch/FLOAT(MAX(ch)) 			; normalize to unity

pctsize=SIZE(pctile)
IF pctsize(0) EQ 0 THEN BEGIN
    vals=LONARR(1)
    pctile=[pctile]
    ENDIF ELSE vals=LONARR(pctsize(N_ELEMENTS(pctsize)-1))
IF KEYWORD_SET(flag) THEN BEGIN
    IF flag EQ 1 $
	THEN FOR k=0,N_ELEMENTS(pctile)-1 DO $
	         vals(k)=hbin*MAX(WHERE(nch LE pctile(k))) $
        ELSE FOR k=0,N_ELEMENTS(pctile)-1 DO $
	         vals(k)=hbin*MIN(WHERE(nch GE pctile(k)))
    ENDIF ELSE BEGIN
    FOR k=0,N_ELEMENTS(pctile)-1 DO $
        IF pctile(k) GE 0.5  $
	    THEN vals(k)=hbin*MAX(WHERE(nch LE pctile(k))) $
            ELSE vals(k)=hbin*MIN(WHERE(nch GE pctile(k)))
    ENDELSE

IF pctsize(0) EQ 0 THEN vals=vals(0)

RETURN,vals

END