website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/d1.f

   1 c===========================================================================
   2 c
   3 c   This file is part of TISEAN
   4 c
   5 c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
   6 c
   7 c   TISEAN is free software; you can redistribute it and/or modify
   8 c   it under the terms of the GNU General Public License as published by
   9 c   the Free Software Foundation; either version 2 of the License, or
  10 c   (at your option) any later version.
  11 c
  12 c   TISEAN is distributed in the hope that it will be useful,
  13 c   but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 c   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 c   GNU General Public License for more details.
  16 c
  17 c   You should have received a copy of the GNU General Public License
  18 c   along with TISEAN; if not, write to the Free Software
  19 c   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20 c
  21 c===========================================================================
  22 c
  23 c   d1 with finite sample correction following Grassberger
  24 c   subroutine for c1
  25 c
  26 c===========================================================================
  27       subroutine d1(nmax,mmax,nxx,y,id,m,ncmin,pr,pln,eln,nmin,kmax)
  28       parameter(im=100,ii=100000000,nx=100000,tiny=1e-20)
  29       dimension y(nxx,mmax),jh(0:im*im),ju(nx),d(nx),jpntr(nx),
  30      .   nlist(nx),nwork(nx)
  31       external rand
  32
  33       if(nmax.gt.nx) stop "d1: make nx larger."
  34       mt=(m-1)/mmax+1
  35       ncomp=nmax-(mt-1)*id
  36       kpr=int(exp(pr)*(ncomp-2*nmin-1))+1
  37       k=int(exp(pln)*(ncomp-2*nmin-1))+1
  38       if(k.gt.kmax) then
  39          ncomp=real(ncomp-2*nmin-1)*real(kmax)/k+2*nmin+1
  40          k=kmax
  41       endif
  42       pln=psi(k)-log(real(ncomp-2*nmin-1))
  43       if(k.eq.kpr) return
  44       write(istderr(),*) 'Mass ', exp(pln),': k=', k, ', N=', ncomp
  45       call rms(nmax,y,sc,sd)
  46       eps=exp(pln/m)*sd
  47       do 10 i=1,nmax-(mt-1)*id
  48  10      ju(i)=i+(mt-1)*id
  49       do 20 i=1,nmax-(mt-1)*id
  50          iperm=min(int(rand(0.0)*nmax-(mt-1)*id)+1,nmax-(mt-1)*id)
  51          ih=ju(i)
  52          ju(i)=ju(iperm)
  53  20      ju(iperm)=ih
  54       iu=ncmin
  55       eln=0
  56  1    call mbase(ncomp+(mt-1)*id,mmax,nxx,y,id,m,jh,jpntr,eps)
  57       iunp=0
  58       do 30 nn=1,iu                                           ! find neighbours
  59          n=ju(nn)
  60          call mneigh(nmax,mmax,nxx,y,n,nmax,id,m,jh,jpntr,eps,
  61      .      nlist,nfound)
  62          nf=0
  63          do 40 ip=1,nfound
  64             np=nlist(ip)
  65             nmd=mod(abs(np-n),ncomp)
  66             if(nmd.le.nmin.or.nmd.ge.ncomp-nmin) goto 40  ! temporal neighbours
  67             nf=nf+1
  68             dis=0
  69             mcount=0
  70             do 50 i=mt-1,0,-1
  71                do 50 is=1,mmax
  72                   mcount=mcount+1
  73                   if(mcount.gt.m) goto 2
  74  50               dis=max(dis,abs(y(n-i*id,is)-y(np-i*id,is)))
  75  2          d(nf)=dis
  76  40         continue
  77          if(nf.lt.k) then
  78             iunp=iunp+1                                   ! mark for next sweep
  79             ju(iunp)=n
  80          else
  81             e=which(nf,d,k,nwork)
  82             eln=eln+log(max(e,tiny))
  83          endif
  84  30      continue
  85       iu=iunp
  86       eps=eps*sqrt(2.)
  87       if(iunp.ne.0) goto 1
  88       eln=eln/(ncmin-(mt-1)*id)
  89       end
  90
  91 c digamma function
  92 c Copyright (C) T. Schreiber (1998)
  93
  94       function psi(i)
  95       dimension p(0:20)
  96       data p/0.,
  97      .  -0.57721566490,  0.42278433509,  0.92278433509,  1.25611766843,
  98      .   1.50611766843,  1.70611766843,  1.87278433509,  2.01564147795,
  99      .   2.14064147795,  2.25175258906,  2.35175258906,  2.44266167997,
 100      .   2.52599501330,  2.60291809023,  2.67434666166,  2.74101332832,
 101      .   2.80351332832,  2.86233685773,  2.91789241329,  2.97052399224/
 102
 103       if(i.le.20) then
 104          psi=p(i)
 105       else
 106          psi=log(real(i))-1/(2.*i)
 107       endif
 108       end
 109
 110
 111
 112
 113
 114
 115
 116