function mkp2correlation,indts,depts,remts,t,filter=filter,refperiod=refperiod,$
  datathresh=datathresh
;
; THIS WORKS WITH REMTS BEING A 2D ARRAY (nseries,ntime) OF MULTIPLE TIMESERIES
; WHOSE INFLUENCE IS TO BE REMOVED.  UNFORTUNATELY THE IDL5.4 p_correlate
; FAILS WITH >1 SERIES TO HOLD CONSTANT, SO I HAVE TO REMOVE THEIR INFLUENCE
; FROM BOTH INDTS AND DEPTS USING MULTIPLE LINEAR REGRESSION AND THEN USE THE
; USUAL correlate FUNCTION ON THE RESIDUALS.
;
; 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'
rsize=size(remts)
case rsize[0] of
  1: begin
    if nt ne n_elements(remts) then message,'indts and remts must be the same size'
    nser=1
    remts=reform(remts,1,nt)
  end
  2: begin
    nser=rsize[1]
    if nt ne rsize[2] then message,'indts and remts must be the same size'
  end
  else: message,'remts must be 1D or 2D'
endcase
;
; 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 multiple series
  remtsh=remts
  remtsl=remts
  for iser = 0 , nser-1 do begin
    ts1=reform(remts[iser,*])
    filter_cru,filter,tsin=ts1,/nan,tshigh=tsh1,tslow=tsl1
    remtsh[iser,*]=tsh1
    remtsl[iser,*]=tsl1
  endfor
  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)
;
remtot=total(remtsr,1)
kl=where(finite(indtsr+deptsr+remtot),nkeep)
;
if nkeep lt dthresh then begin
  oner=!values.f_nan
endif else begin
  ; IDL5.4 p_correlate FAILS WHEN nser>1
  ;oner=p_correlate(indtsr(kl),deptsr(kl),reform(remtsr(kl),1,nkeep))
  ; Instead use MLR to remove influences from BOTH ind and dep timeseries
  dummy=regress(remtsr[*,kl],indtsr[kl],yfit=ifit)
  dummy=regress(remtsr[*,kl],deptsr[kl],yfit=dfit)
  oner=correlate(indtsr[kl]-ifit[*],deptsr[kl]-dfit[*])
endelse
;
if dofilt eq 1 then begin
  ;
  if nkeep lt dthresh then begin
    oner=replicate(!values.f_nan,3)
  endif else begin
    dummy=regress(remtsh[*,kl],indtsh[kl],yfit=ifit)
    dummy=regress(remtsh[*,kl],deptsh[kl],yfit=dfit)
    oner=[oner,correlate(indtsh[kl]-ifit[*],deptsh[kl]-dfit[*])]
    dummy=regress(remtsl[*,kl],indtsl[kl],yfit=ifit)
    dummy=regress(remtsl[*,kl],deptsl[kl],yfit=dfit)
    oner=[oner,correlate(indtsl[kl]-ifit[*],deptsl[kl]-dfit[*])]
    ;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