1 c===========================================================================
3 c This file is part of TISEAN
5 c Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
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.
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.
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
21 c===========================================================================
23 c d1 with finite sample correction following Grassberger
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),
33 if(nmax.gt.nx) stop "d1: make nx larger."
36 kpr=int(exp(pr)*(ncomp-2*nmin-1))+1
37 k=int(exp(pln)*(ncomp-2*nmin-1))+1
39 ncomp=real(ncomp-2*nmin-1)*real(kmax)/k+2*nmin+1
42 pln=psi(k)-log(real(ncomp-2*nmin-1))
44 write(istderr(),*) 'Mass ', exp(pln),': k=', k, ', N=', ncomp
45 call rms(nmax,y,sc,sd)
47 do 10 i=1,nmax-(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)
56 1 call mbase(ncomp+(mt-1)*id,mmax,nxx,y,id,m,jh,jpntr,eps)
58 do 30 nn=1,iu ! find neighbours
60 call mneigh(nmax,mmax,nxx,y,n,nmax,id,m,jh,jpntr,eps,
65 nmd=mod(abs(np-n),ncomp)
66 if(nmd.le.nmin.or.nmd.ge.ncomp-nmin) goto 40 ! temporal neighbours
73 if(mcount.gt.m) goto 2
74 50 dis=max(dis,abs(y(n-i*id,is)-y(np-i*id,is)))
78 iunp=iunp+1 ! mark for next sweep
82 eln=eln+log(max(e,tiny))
88 eln=eln/(ncmin-(mt-1)*id)
92 c Copyright (C) T. Schreiber (1998)
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/
106 psi=log(real(i))-1/(2.*i)