function mkpcorrelation,indts,depts,remts,t,filter=filter,refperiod=refperiod,$
  datathresh=datathresh
;
; DO NOT EXTEND THIS TO WORK WITH MULTIPLE remts TIMESERIES, BECAUSE IDL5.4's
; p_correlate FUNCTION IS IN ERROR WHEN MULTIPLE SERIES ARE USED!
;
; Given three timeseries, indts, depts and remts, possibly containing
; missing data given by NaN, this function returns the correlation between
; the dependent timeseries (depts) and the independent timeseries (indts),
; having first removed the influence of remts on indts, providing
; there are at least datathresh (or 5 if datathresh is not specified) triples
; of non-missing values.  If datathresh is specified, but is <= 1,
; then it is a fraction of the length of the input timeseries.  The
; correlation is computed over the entire timeseries (which must be equal
; length), unless limited by a given refperiod (assumed to be in elements
; numbers, unless a time variable is also given in t).  If filter is
; specified, then a 3-element vector is returned contained the correlations
; between the raw timeseries, the high-pass timeseries and the low-pass
; timeseries, with the filtering done via a gaussian-weighted filter with
; half-width given by the value of filter.
;
;-------------------------------------------------------------------
;
; Get input series sizes
;
nt=n_elements(indts)
if nt ne n_elements(depts) then message,'indts and depts must be the same size'
if nt ne n_elements(remts) then message,'indts and remts must be the same size'
;
; Create or check time values
;
if n_params() eq 2 then begin
  t=findgen(nt)
endif else begin
  if n_elements(t) ne nt then message,'t must be same size as timeseries'
endelse
;
; Create or check reference period
;
case n_elements(refperiod) of
  0: begin
    startel=0
    endel=nt-1
    end
  2: begin
    startel=where(t eq refperiod(0),n)
    if n ne 1 then message,'Cannot find the ref start time in the t vector'
    endel=where(t eq refperiod(1),n)
    if n ne 1 then message,'Cannot find the ref end time in the t vector'
    startel=startel(0)
    endel=endel(0)
    end
  else: message,'refperiod must be a 2-element [start,end] vector'
endcase
;
; Compute minimum data requirements
;
if n_elements(datathresh) eq 0 then begin
  dthresh=5.
endif else begin
  if datathresh le 1. then begin
    dthresh=datathresh*(endel-startel+1)
    dthresh=dthresh > 5.
  endif else begin
    dthresh=datathresh
    dthresh=dthresh < (endel-startel+1)
    dthresh=dthresh > 5.
  endelse
endelse
;
dofilt=0
if n_elements(filter) ne 0 then begin
  dofilt=1
  filter_cru,filter,tsin=indts,/nan,tshigh=indtsh,tslow=indtsl
  filter_cru,filter,tsin=depts,/nan,tshigh=deptsh,tslow=deptsl
  filter_cru,filter,tsin=remts,/nan,tshigh=remtsh,tslow=remtsl
  indtsh=indtsh(startel:endel)
  indtsl=indtsl(startel:endel)
  deptsh=deptsh(startel:endel)
  deptsl=deptsl(startel:endel)
  remtsh=remtsh(startel:endel)
  remtsl=remtsl(startel:endel)
endif
indtsr=indts(startel:endel)
deptsr=depts(startel:endel)
remtsr=remts(startel:endel)
;
kl=where(finite(indtsr+deptsr+remtsr),nkeep)
;
if nkeep lt dthresh then oner=!values.f_nan $
  else oner=p_correlate(indtsr(kl),deptsr(kl),reform(remtsr(kl),1,nkeep))
;
if dofilt eq 1 then begin
  ;
  if nkeep lt dthresh then begin
    oner=replicate(!values.f_nan,3)
  endif else begin
    oner=[oner,p_correlate(indtsh(kl),deptsh(kl),reform(remtsh(kl),1,nkeep))]
    oner=[oner,p_correlate(indtsl(kl),deptsl(kl),reform(remtsl(kl),1,nkeep))]
  endelse
  ;
endif
;
return,oner
;
end