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