website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/c1.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   information dimension, fixed mass
  23 c   see  H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, Cambridge
  24 c      University Press (1997)
  25 c   author T. Schreiber (1999)
  26 c===========================================================================
  27
  28       parameter(nx=100000,mx=10)
  29       dimension x(nx,mx), icol(mx)
  30       character*72 file, fout
  31       data kmax/100/, res/2./
  32       data iverb/1/
  33       external rand
  34
  35       call whatido("fixed mass approach to d1 estimation",iverb)
  36       id=imust("d")
  37       mfrom=imust("m")
  38       mto=imust("M")
  39       ntmin=imust("t")
  40       ncmin=imust("n")
  41       res=fcan("#",res)
  42       r=rand(sqrt(abs(fcan("I",0.0))))
  43       kmax=ican("K",kmax)
  44       resl=log(2.)/res
  45       nmax=ican("l",nx)
  46       nexcl=ican("x",0)
  47       call columns(mc,mx,icol)
  48       mcmax=max(1,mc)
  49       isout=igetout(fout,iverb)
  50       if(fout.eq." ") isout=1
  51
  52       call nthstring(1,file)
  53       call xreadfile(nmax,mcmax,nx,x,nexcl,icol,file,iverb)
  54       if(file.eq."-") file="stdin"
  55       if(isout.eq.1) call addsuff(fout,file,"_c1")
  56       call outfile(fout,iunit,iverb)
  57       do 10 m=mfrom,mto
  58          write(iunit,'(4h#m= ,i5)') m
  59          pr=0.
  60          do 20 pl=log(1./(nmax-(m-1)*id)),0.,resl
  61             pln=pl
  62             call d1(nmax,mcmax,nx,x,id,m,ncmin,pr,pln,rln,ntmin,kmax)
  63             if(pln.eq.pr) goto 20
  64             it=it+1
  65             pr=pln
  66             write(iunit,*)   exp(rln), exp(pln)
  67  20         continue
  68          write(iunit,'()')
  69  10      write(iunit,'()')
  70       end
  71
  72       subroutine usage()
  73 c usage message
  74
  75       call whatineed(
  76      .   "-d# -m# -M# -t# -n# "//
  77      .   "[-## -K# -o outfile -I# -l# -x# -c#,# -V# -h] file")
  78       call popt("d","delay")
  79       call popt("m","minimal total embedding dimension")
  80       call popt("M","maximal total embedding dimension")
  81       call popt("t","minimal time separation")
  82       call popt("n","minimal number of center points")
  83       call popt("#","resolution, values per octave (2)")
  84       call popt("K","maximal number of neighbours (100)")
  85       call popt("I","seed for random numbers")
  86       call popt("l","number of values to be read (all)")
  87       call popt("x","number of values to be skipped (0)")
  88       call popt("c","column(s) to be read (1 or file,#)")
  89       call pout("file_c1")
  90       call pall()
  91       stop
  92       end