function fclazy(nmax,y,id,m,s,ifo,eps0)
      dimension y(nmax),s(m)

      eps=eps0
      nfound=0
 1    av=0
      do 10 n=(m-1)*id+1,nmax-ifo
        do 20 i=1,m
 20       if(abs(y(n-(m-i)*id)-s(i)).ge.eps) goto 10
        nfound=nfound+1                                ! if got here: neighbour
        av=av+y(n+ifo)                      ! average over future of neighbours
 10     continue
      eps=eps*1.2
      if(nfound.eq.0) goto 1                         ! try larger neighbourhood
      fclazy=av/nfound
      end

      subroutine base(nmax,y,id,m,jh,jpntr,eps)
      parameter(im=100,ii=100000000) 
      dimension y(nmax),jh(0:im*im),jpntr(nmax)

      do 10 i=0,im*im
 10     jh(i)=0
      do 20 n=(m-1)*id+1,nmax                                  ! make histogram
        i=im*mod(int(y(n)/eps)+ii,im)+mod(int(y(n-(m-1)*id)/eps)+ii,im)
 20     jh(i)=jh(i)+1
      do 30 i=1,im*im                                           ! accumulate it
 30     jh(i)=jh(i)+jh(i-1)
      do 40 n=(m-1)*id+1,nmax                           ! fill list of pointers
        i=im*mod(int(y(n)/eps)+ii,im)+mod(int(y(n-(m-1)*id)/eps)+ii,im)
        jpntr(jh(i))=n
 40     jh(i)=jh(i)-1
      end

      subroutine neigh(nmax,y,n,nlast,id,m,jh,jpntr,eps,nlist,nfound)
      parameter(im=100,ii=100000000) 
      dimension y(nmax),jh(0:im*im),jpntr(nmax),nlist(nmax)

      nfound=0
      jj=int(y(n)/eps)
      kk=int(y(n-(m-1)*id)/eps)
      do 10 j=jj-1,jj+1                               ! scan neighbouring boxes
         do 20 k=kk-1,kk+1
            jk=im*mod(j+ii,im)+mod(k+ii,im)
            do 30 ip=jh(jk+1),jh(jk)+1,-1               ! this is in time order
               np=jpntr(ip)
               if(np.gt.nlast) goto 20
               do 40 i=0,m-1
 40               if(abs(y(n-i*id)-y(np-i*id)).ge.eps) goto 30
               nfound=nfound+1
               nlist(nfound)=np                       ! make list of neighbours
 30            continue
 20         continue
 10      continue
      end

      subroutine nrlazy(nmax,y,yc,m,eps)
      parameter(im=100,ii=100000000,nx=50000) 
      dimension y(nmax),yc(nmax),jh(0:im*im),jpntr(nx),nlist(nx)

      if(nmax.gt.nx) pause "Make nx larger."
      call base(nmax,y,1,m,jh,jpntr,eps)
      do 10 n=1,nmax
 10     yc(n)=y(n)   
      do 20 n=m,nmax           
        call neigh(nmax,y,n,nmax,1,m,jh,jpntr,eps,nlist,nfound)
        av=0
        do 30 nn=1,nfound            
 30       av=av+y(nlist(nn)-(m-1)/2)                ! average middle coordinate
 20     yc(n-(m-1)/2)=av/nfound
      end

      subroutine lyap(nmax,y,id,m,eps,ifu,s,nmin,nfmin,ncmin)
      parameter (im=100,ii=100000000,nx=50000,ifum=200,tiny=1e-20)
      dimension y(nmax),s(0:ifu),sh(0:ifum),jh(0:im*im),
     .  jpntr(nx),nlist(nx)

      if(nmax.gt.nx) pause "Make nx larger."
      call base(nmax-ifu,y,id,m,jh,jpntr,eps)
      do 10 i=0,ifu
 10     s(i)=0
      nc=0
      do 20 n=(m-1)*id+1, nmax-ifu                           ! reference points
        call neigh(nmax-ifu,y,n,nmax,id,m,jh,jpntr,eps,nlist,nfound)
        do 30 i=0,ifu
 30       sh(i)=0
        nf=0
        do 40 nn=1,nfound
          np=nlist(nn)
          if(abs(n-np).le.nmin) goto 40
          nf=nf+1
          do 50 i=0,ifu                                     ! average distances
 50         sh(i)=sh(i)+abs(y(n+i)-y(np+i))                  
 40       continue
        if(nf.ge.nfmin) then              ! enough neighbours closer $\epsilon$
          nc=nc+1
          do 60 i=0,ifu                     ! average log of averaged distances
 60         s(i)=s(i)+log(max(sh(i)/nf,tiny))  
        endif
 20     if(nc.eq.ncmin) goto 1
 1    if(nc.eq.0) return                     ! no points with enough neighbours
      do 70 i=0,ifu
 70     s(i)=s(i)/nc
      print*, 'tried ', n,' reference points, found ', nc   !c
      end

      subroutine correl(nmax,y,eps,id,m,c,nmin,ncmin,ipmin)
      parameter(im=100,ii=100000000,nx=1000000,mm=15)
      dimension y(nmax),jh(0:im*im),ipairs(mm),c(m),jpntr(nx),nlist(nx)

      if(nmax.gt.nx.or.m.gt.mm) pause "Make mm/nx larger."
      do 10 i=1,m-1
 10     ipairs(i)=0
      call base(nmax,y,id,2,jh,jpntr,eps)
      do 20 n=nmin+(m-1)*id+1,nmax
        call neigh(nmax,y,n,n-nmin,id,2,jh,jpntr,eps,nlist,nfound)
        ipairs(1)=ipairs(1)+nfound
        do 30 nn=1,nfound                    ! all neighbours in two dimensions
          np=nlist(nn)
          if(np.lt.(m-1)*id+1) goto 30
          do 40 i=2,m-1
            if(abs(y(n-i*id)-y(np-i*id)).ge.eps) goto 30
 40         ipairs(i)=ipairs(i)+1       ! neighbours in $3,\ldots,m$ dimensions
 30       continue
 20     if(n-nmin-(m-1)*id.ge.ncmin.and.ipairs(m-1).ge.ipmin) goto 1
 1    s=real(n-nmin-(m-1)*id+1)*real(n-nmin-(m-1)*id)/2.        ! normalisation
      do 70 i=1,m-1
 70     c(i+1)=ipairs(i)/s
      end

      subroutine stplot(nmax,y,id,m,epsmax,stp,mfrac,mdt,idt)
      parameter(meps=1000)
      dimension y(nmax),stp(mfrac,mdt),ihist(meps)

      do 10 it=1,mdt
        do 20 ieps=1,meps
 20       ihist(ieps)=0
        print*, it*idt           !c
        do 30 n=it*idt+(m-1)*id+1,nmax
          dis=0                            ! compute distance in $m$ dimensions
          do 40 me=0,m-1
 40         dis=max(dis,abs(y(n-me*id)-y(n-me*id-it*idt)))
          ih=min(int(meps*dis/epsmax)+1,meps)
 30       ihist(ih)=ihist(ih)+1
        do 10 ifrac=1,mfrac
          need=(nmax-it*idt-(m-1)*id)*ifrac/real(mfrac)
          is=0
          do 50 ieps=1,meps
            is=is+ihist(ieps)
 50         if(is.ge.need) goto 1
 1        stp(ifrac,it)=ieps*epsmax/meps
 10     continue
      end

      subroutine clean(nmax,y,yc,m,kmin,nq,d)
      parameter(im=100,ii=100000000,nx=50000,mm=15,small=0.0001) 
      dimension y(nmax),yc(nmax),r(mm),ju(nx),c(mm,mm),cm(mm),
     .  jh(0:im*im),jpntr(nx),nlist(nx)
      dimension zcm(mm,nx)                                                ! \TR

      if(nmax.gt.nx.or.m.gt.mm) pause "Make mm/nx larger."
      sr=2*small+m-2                                        ! ${\rm tr}(1/r)=1$
      do 10 i=1,m
        r(i)=sr
 10     if(i.eq.m.or.i.eq.1) r(i)=sr/small
      do 20 i=1,nmax
 20     yc(i)=y(i)
      do 30 istep=1,2                                                     ! \TR
        eps=d
        iu=nmax-m+1
        do 40 i=1,iu
 40        ju(i)=i+m-1
 1      call base(nmax,y,1,m,jh,jpntr,eps)
        iunp=0
        do 50 nn=1,iu                                         ! find neighbours
          n=ju(nn)
          call neigh(nmax,y,n,nmax,1,m,jh,jpntr,eps,nlist,nfound)
          if(nfound.lt.kmin) then                 ! not enough neighbours found
            iunp=iunp+1                                   ! mark for next sweep
            ju(iunp)=n
          else                                        ! fine: enough neighbours
            do 90 i=1,m                                 ! centre of mass vector
              s=0
              do 100 np=1,nfound
 100            s=s+y(nlist(np)-m+i)
 90           cm(i)=s/nfound
            if(istep.eq.1) then               ! just store centre of mass ! \TR
              do 110 i=1,m                                                ! \TR
 110            zcm(i,n)=cm(i)                                            ! \TR
            else                                                          ! \TR
              do 120 i=1,m              ! corrected centre of mass vector ! \TR
                s=0                                                       ! \TR
                do 130 np=1,nfound                                        ! \TR
 130              s=s+zcm(i,nlist(np))                                    ! \TR
 120            cm(i)=2*cm(i)-s/nfound                                    ! \TR
              do 140 i=1,m                          ! compute covariance matrix
                do 140 j=i,m
                  s=0
                  do 150 np=1,nfound
                    jm=nlist(np)-m
 150                s=s+(y(jm+i)-cm(i))*(y(jm+j)-cm(j))
                  c(i,j)=r(i)*r(j)*s/nfound
 140              c(j,i)=c(i,j)
              call eigen(c,m)                  ! find eigenvectors (increasing)
              do 160 i=1,m
                s=0
                do 170 iq=1,nq
                  do 170 j=1,m
 170                s=s+(y(n-m+j)-cm(j))*c(i,iq)*c(j,iq)*r(j)
 160            yc(n-m+i)=yc(n-m+i)-s/r(i)/r(i)
            endif
          endif
 50       continue
        iu=iunp
        print*, "With ", eps, iunp, " uncorrected" !c
        eps=eps*sqrt(2.)
 30     if(iunp.ne.0) goto 1
      end

      subroutine d1(nmax,y,id,m,ncmin,pln,eln,nmin,kmax)
      parameter(im=100,ii=100000000,nx=50000,tiny=1e-20) 
      dimension y(nmax),jh(0:im*im),ju(nx),d(nx),jpntr(nx),nlist(nx)

      if(nmax.gt.nx) pause "Make nx larger."
      ncomp=nmax-(m-1)*id
      k=int(exp(pln)*(ncomp-2*nmin-1))+1
      if(k.gt.kmax) then
        ncomp=real(ncomp-2*nmin-1)*real(kmax)/k+2*nmin+1
        k=kmax
      endif         
      pln=psi(k)-log(real(ncomp-2*nmin-1))
      print*, 'Mass ', exp(pln),': k=', k, ', N=', ncomp     !c
      eps=exp(pln/m)*sigma(nmax,y)
      iu=ncmin-(m-1)*id
      do 10 i=1,iu
 10     ju(i)=i+(m-1)*id
      eln=0
 1    call base(ncomp+(m-1)*id,y,id,m,jh,jpntr,eps)
      iunp=0
      do 20 nn=1,iu                                           ! find neighbours
        n=ju(nn)
        call neigh(nmax,y,n,nmax,id,m,jh,jpntr,eps,nlist,nfound)
        nf=0
        do 30 ip=1,nfound
          np=nlist(ip)
          nmd=mod(abs(np-n),ncomp)
          if(nmd.le.nmin.or.nmd.ge.ncomp-nmin) goto 30    ! temporal neighbours
          nf=nf+1
          dis=0
          do 40 i=0,m-1
 40         dis=max(dis,abs(y(n-i*id)-y(np-i*id)))
          d(nf)=dis
 30       continue
        if(nf.lt.k) then
          iunp=iunp+1                                     ! mark for next sweep
          ju(iunp)=n
        else
          e=which(k,nf,d)
          eln=eln+log(max(e,tiny))
        endif
 20     continue
      iu=iunp
      print*, eps, iunp  !c
      eps=eps*sqrt(2.)
      if(iunp.ne.0) goto 1
      eln=eln/(ncmin-(m-1)*id)
      end