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 correlation integral c2, no fast neighbour search
24 c complete direct neighbour search
25 c author T. Schreiber (1998)
26 c modified H. Kantz, Feb. 2007
27 c===========================================================================
29 parameter(nx=1000000,me=30,meps=800)
30 dimension x(nx), c(0:meps,me)
31 c j=a*(log(d)-log(xmax-xmin)) d=xmax-xmin -> j=0
32 c a=-rs/log(2.) d=s*2**(-j/res)
33 character*72 file, fout
37 call whatido("correlation sum, complete naive neighbour search,
38 .univariate data only"
40 call whatido("univariate data only",iverb)
50 isout=igetout(fout,iverb)
51 if(fout.eq." ") isout=1
53 do 10 ifi=1,nstrings()
54 call nthstring(ifi,file)
56 call readfile(nmax,x,nexcl,jcol,file,iverb)
57 if(file.eq."-") file="stdin"
58 if(isout.eq.1) call addsuff(fout,file,"_c2")
59 call minmax(nmax,x,xmin,xmax)
65 call d2naive(nmax,x,id,mmin,mmax,c,meps,log(sc),a,ntmin,ntmax)
66 call outfile(fout,iunit,iverb)
68 write(iunit,'(4h#m= ,i5)') m
70 40 c(j,m)=c(j,m)+c(j+1,m)
72 if(c(j,m).eq.0.) goto 1
73 50 write(iunit,*) sc*2**(-j/res), c(j,m)/c(0,m)
84 . "-d# -M# -t# [-m# -##"//
85 . " -o outfile -l# -x# -c# -V# -h] file(s)")
86 call popt("d","delay")
87 call popt("M","maximal embedding dimension")
88 call popt("t","minimal time separation")
89 call popt("m","minimal embedding dimension [1]")
90 call popt("#","resolution, values per octave [2]")
91 c call popt("T","for Guido")
92 call popt("l","number of values to be read [all]")
93 call popt("x","number of values to be skipped [0]")
94 call popt("c","column to be read [1 or file,#]")
100 subroutine d2naive(nmax,x,id,mmin,mmax,c,meps,scl,a,ntmin,ntmax)
101 parameter(nx=1000000,tiny=1e-30)
102 dimension x(nmax),c(0:meps,mmax),d(nx)
104 if(nmax.gt.nx) stop "d2naive: make nx larger."
105 nlast=min(nmax-(mmax-1)*id-1,ntmax)
106 do 10 ndt=ntmin,nlast
108 20 d(n)=max(abs(x(n)-x(n-ndt)),tiny)
109 do 10 n=ndt+1+(mmax-1)*id,nmax
112 30 dmax=max(dmax,d(n-(m-1)*id))
113 j=int(a*(log(dmax)-scl))
115 if(d(n-(m-1)*id).gt.dmax) then
117 j=int(a*(log(dmax)-scl))