pro distrib_trendmean,tsin,timey,seglen=seglen,$
  binmean=binmean,bmhalf=bmhalf,bintren=bintren,bthalf=bthalf,$
  tren95=tren95,mean95=mean95,freqmean=freqmean,freqtren=freqtren
;
; Given an input time series, this function computes the distributions
; and 95% cutoffs of trends and means (from overlapping segments) with
; the given segment length.  The centres of the bins are given (should
; be equally-spaced, and end bins are open), as is the half-width of
; the bins.
;
nbin=n_elements(binmean)
if n_elements(bintren) ne nbin then message,'Ooops!'
freqmean=fltarr(nbin)
freqtren=fltarr(nbin)
;
nyr=n_elements(tsin)
ntren=nyr-seglen+1
obstren=fltarr(ntren)*!values.f_nan
obsmean=fltarr(ntren)*!values.f_nan
;
; Process each overlapping segment separately, updating the histograms
; of trends and means as we go
for k = 0 , ntren-1 do begin
  xval=timey(k:k+seglen-1)
  yval=tsin(k:k+seglen-1)
  kl=where(finite(yval),nkeep)
  if nkeep gt 0.66*seglen then begin     ; need at least 2/3rds data
    xval=xval(kl)
    yval=yval(kl)
    acoeff=linfit(xval,yval)
    obstren(k)=acoeff(1)*float(seglen)
    obsmean(k)=total(yval)/float(nkeep)
    ; Identify which bins these values fall in
    kkk=where(obsmean(k) ge (binmean-bmhalf),nkkk)
    if nkkk ge 1 then kkk=kkk(nkkk-1) else kkk=0
    freqmean(kkk)=freqmean(kkk)+1
    kkk=where(obstren(k) ge (bintren-bthalf),nkkk)
    if nkkk ge 1 then kkk=kkk(nkkk-1) else kkk=0
    freqtren(kkk)=freqtren(kkk)+1
  endif
endfor
kl=where(finite(obstren),ntren)
obstren=obstren(kl)
kl=where(finite(obsmean))
obsmean=obsmean(kl)
print,'No. of means or trends analysed:',ntren
;
; Compute the 95% threshold (easy cos there's no missing data)
sortobs=obstren(sort(obstren))
i95=ntren*0.95
tren95=sortobs(i95)
sortobs=obsmean(sort(obsmean))
i95=ntren*0.95
mean95=sortobs(i95)
;
end