website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/ar-run.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   ar-run.f
  23 c   iterate AR model, e.g. as fitted by ar-model (Dresden)
  24 c   author T. Schreiber (1999)
  25 c===========================================================================
  26
  27       parameter(npmax=100)
  28       character*72 file, fout, fline
  29       dimension x(-npmax:npmax), a(npmax)
  30       external rand
  31       data np/npmax/, ntrans/10000/, iuni/0/
  32       data iverb/1/
  33
  34       call whatido("iterate AR model, e.g. as fitted by ar-model",iverb)
  35       np=ican("p",np)
  36       if(np.gt.npmax) stop "ar-run: make npmax larger."
  37       nmax=imust('l')
  38       ntrans=ican("x",ntrans)
  39       if(lopt("u",1).eq.1) iuni=1
  40       r=rand(sqrt(abs(fcan("I",0.0))))
  41       isout=igetout(fout,iverb)
  42
  43       do 10 n=1,npmax
  44          x(-n)=0.
  45  10      x(n)=0.
  46       call nthstring(1,file)
  47       call infile(file,iunit,iverb)
  48       read(iunit,'(a)') fline
  49       if(fline(1:1).eq."#") then
  50          read(fline(18:72),'(f20.0)',err=999) var
  51          do 20 j=1,np
  52             read(iunit,'(a1,f20.0)',err=999) fline(1:1), a(j)
  53  20         if(fline(1:1).ne."#") goto 1
  54       else
  55          read(fline(1:72),'(f20.0)',err=999) var
  56          do 30 j=1,np
  57  30         read(iunit,'(f20.0)',err=999,end=1) a(j)
  58       endif
  59  1    np=j-1
  60       if(iv_echo(iverb).eq.1) then
  61          write(istderr(),*) 'coefficients:      ', (a(i),i=1,np)
  62          write(istderr(),*) 'driving amplitude: ', var
  63       endif
  64       if(isout.eq.1) fout="ar.dat"
  65       call outfile(fout,iunit,iverb)
  66       n=-ntrans
  67  2    n=n+1
  68       nn=mod(n+ntrans,np)+1
  69       xx=rgauss(0.0,var)
  70       do 40 j=1,np
  71  40      xx=xx+a(j)*x(nn-j)
  72       x(nn)=xx
  73       x(nn-np)=xx
  74       if(n.lt.1) goto 2
  75       write(iunit,*) xx
  76       if(nmax.eq.0.or.n.lt.nmax) goto 2
  77       stop
  78
  79  999  write(istderr(),'(a)') "wrong input format! try:"
  80       write(istderr(),'(a)') "(rms of increments)"
  81       write(istderr(),'(a)') "a(1)"
  82       write(istderr(),'(a)') "a(2)"
  83       write(istderr(),'(a)') "..."
  84       end
  85
  86       subroutine usage()
  87 c usage message
  88
  89       call whatineed(
  90      .   "-l# [-p# -I# -o outfile -x# -V# -h] file")
  91       call popt("l","number of iterations (l=0: infinite)")
  92       call popt("p","order of AR-model (default determined from input)")
  93       call popt("I","seed for random numbers")
  94       call popt("x","number of transients discarded (10000)")
  95       call pout("ar.dat")
  96       call pall()
  97       stop
  98       end
  99