;Find the value which marks the bottom q% of data
;
;Input
;  data - Data to find the quantile of. Should be an array, may be unsorted
;  q - Fraction of data. May be a scalar or array, but all values should be
;      between 0 and 1, inclusive.
;  /nan - don't include NaN values. Otherwise, they fall where the sort puts them.
;Output
;  Value for each quantile in q. Corresponding results in q and return have
;  the same index
;
;Simple Example:
;  Take a simple index array of data.
;    > a=findgen(1000)
;  Naturally, the top 10% of this data is the data above 900, so we want
;  quantile() to tell us that
;    > print,quantile(a,0.9)
;       900.000
;
;More complex example:
;  What are the levels on a sine wave that 5% of the points are below and 5% are above?
;   > a=sin(2*!pi*findgen(10000)/10000.0)
;   > print,quantile(a,[0.05,0.95])
;    -0.987688     0.987688
function quantile,data_,q,nan=nan
  if keyword_set(nan) then begin
    w=where(finite(data_),count)
    if count eq 0 then return,!values.f_nan
    data=data_[w]
  end else begin
    data=data_
  end
  s=sort(data)
  qq=q*n_elements(data)
  return,data[s[qq]]
end