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