website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/c2t.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   c2t.f
  23 c   Takens' estimator from c2
  24 c   author T. Schreiber (1998)
  25 c===========================================================================
  26
  27       parameter(meps=1000)
  28       dimension e(meps), c(meps), lw(meps)
  29       double precision a, b
  30       character*72 file, fout, aline
  31       data iverb/1/
  32
  33       call whatido("Takens' estimator from correlation sum data",iverb)
  34       isout=igetout(fout,iverb)
  35       if(nstrings().eq.0) call usage()
  36       call nthstring(1,file)
  37       call infile(file,iunit,iverb)
  38       if(isout.eq.1) call addsuff(fout,file,"_t")
  39       call outfile(fout,iunit2,iverb)
  40  1    read(iunit,'(a)',end=999) aline
  41  4    if(aline(1:1).ne."#") goto 1
  42       if(aline(1:1).eq."#")
  43      .   read(aline(index(aline,"m=")+2:72),'(i20)',err=1) m
  44       me=0
  45  2    read(iunit,'(a)') aline
  46       if(aline(1:72).eq." ") goto 3
  47       read(aline,*,err=999,end=999) ee, cc
  48       if(cc.le.0.) goto 3
  49       me=me+1
  50       e(me)=log(ee)
  51       c(me)=log(cc)
  52       goto 2
  53  3    write(iunit2,'(4h#m= ,i5)') m
  54       call indexx(me,e,lw)
  55       call index2sort(me,e,lw)
  56       call index2sort(me,c,lw)
  57       cint=0
  58       do 10 i=2,me
  59          b=(e(i)*c(i-1)-e(i-1)*c(i))/(e(i)-e(i-1))
  60          a=(c(i)-c(i-1))/(e(i)-e(i-1))
  61          if(a.ne.0) then
  62             cint=cint+(exp(b)/a)*(exp(a*e(i))-exp(a*e(i-1)))
  63          else
  64             cint=cint+exp(b)*(e(i)-e(i-1))
  65          endif
  66  10      write(iunit2,*) exp(e(i)), exp(c(i))/cint
  67       write(iunit2,'()')
  68       write(iunit2,'()')
  69       goto 4
  70  999  stop
  71       end
  72
  73       subroutine usage()
  74 c usage message
  75
  76       call whatineed(
  77      .   "[-o outfile -V# -h] file")
  78       call pout("file_t")
  79       call pall()
  80       stop
  81       end