pro fill_t, reihe, reihe1, neureihe, z

;+
; NAME:
;       FILL_T
; PURPOSE:
;       Filling invalid points outside attractors. Should be used to the output of 'fil_t_att'
; EXPLANATION:
;       readme.txt
;
;
; CALLING SEQUENCE:
;       fill_t, reihe, reihe1, neureihe [,z]
; INPUTS:
;       reihe: Detrended level 1 time series to be filled ((n) vector)
;       riehe1: Detrended level 1 time series used for filling ((n) vector)
; OUTPUTS:
;       neureihe: Filled time series
; PROCEDURES USED:
;
;       Procedures:  WO
;
; MODIFICATION HISTORY:
;       Package 'SPM level 1.5' written between 2000 and 2003, Richard Wachter, PMOD/WRC
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

 neureihe=reihe
 l=n_elements(reihe)
 gblstd=stddev(reihe, /nan)
 gblstd1=stddev(reihe1, /nan)

 wo, reihe, nullen

 lg=(size(nullen))[2]
 stv=replicate(!values.f_nan, lg)
 sth=replicate(!values.f_nan, lg)


 for i=0, lg-1 do begin

  as=2880*2
  wdh: as=as/2

  vh=reihe[(nullen[0,i]-as)>0:nullen[0,i]-1>0]
  hh=reihe[nullen[1,i]+1:(nullen[1,i]+as)<l-1]

  fak=replicate(!values.f_nan, 2)
  z=0
  if n_elements(where(finite(hh))) gt as/2 then begin

    z=z+1
    stv[i]=stddev(hh, /nan)
    fak[0]=stv[i]

  endif


  if n_elements(where(finite(vh))) gt as/2 then begin

   z=z+1
   sth[i]=stddev(vh, /nan)
   fak[1]=sth[i]

  endif

  if z gt 0 then begin


   if abs(fak[1]-gblstd)/gblstd gt 0.18 then fak[1]=!values.f_nan
   if abs(fak[0]-gblstd)/gblstd gt 0.18  then fak[0]=!values.f_nan

  endif


  faktor=mean(fak, /nan)

  if not finite(faktor) then  begin

    if as gt 180 then goto, wdh else print, i, ' problem', nullen[*,i]

  endif else begin

   la=min([long((nullen[1,i]-nullen[0,i]+1)/2.)+1L, 30])

   ug=max([0,nullen[0,i]-la])
   og=min([l-1L,nullen[1,i]+la])

   ind=findgen(la+1)/(la+1)

   seq1=reihe1[ug:og]





   if n_elements(where(finite(seq1))) gt 180 then faks1=stddev(seq1, /nan) else faks1=gblstd1


   nr=reihe1*faktor/faks1

   if ug le nullen[0,i]-1 and (where(abs(nr[ug:nullen[0,i]-1]) lt 0.0))[0] eq -1 $
     then neureihe[ug:nullen[0,i]-1]=ind*nr[ug:nullen[0,i]-1] + (1.-ind)*reihe[ug:nullen[0,i]-1]

   neureihe[nullen[0,i]:nullen[1,i]]=nr[nullen[0,i]:nullen[1,i]]

   if og ge nullen[1,i]+1 and (where(abs(nr[nullen[1,i]+1:og]) lt 0.0))[0] eq -1 $
     then neureihe[nullen[1,i]+1:og]=ind*reihe[nullen[1,i]+1:og] + (1. - ind)*nr[nullen[1,i]+1:og]

  endelse

 endfor

end