binaries/src/disembl/Tisean_3.0.1/source_f/endtoend.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   endtoend.f
  23 c   Determine end-to-end mismatch before making surrogate data
  24 c   author T. Schreiber (1999)
  25 c===========================================================================
  26
  27       parameter(nx=100000,mx=20)
  28       dimension x(nx,mx), icol(mx)
  29       character*72 file, fout
  30       data iverb/15/
  31
  32       call whatido("Determine end-to-end mismatch",iverb)
  33       nmax=ican("l",nx)
  34       nexcl=ican("x",0)
  35       wjump=fcan("j",0.5)
  36       mcmax=ican("m",0)
  37       call columns(mc,mx,icol)
  38       if(mcmax.eq.0) mcmax=max(1,mc)
  39       isout=igetout(fout,iverb)
  40
  41       call nthstring(1,file)
  42       call xreadfile(nmax,mcmax,nx,x,nexcl,icol,file,iverb)
  43       if(file.eq."-") file="stdin"
  44       if(isout.eq.1) call addsuff(fout,file,"_end")
  45       call outfile(fout,iunit,iverb)
  46       nmaxp=nmax
  47       etotm=mcmax
  48  1    nmaxp=nless(nmaxp)
  49       call jump(nmax,nmaxp,nx,x,mcmax,wjump,etot,ejump,eslip,njump)
  50       if(etot.lt.etotm) then
  51          etotm=etot
  52          write(iunit,'(a,i7,a,i7,a,f5.1,a)')
  53      .      "length:", nmaxp,
  54      .      "  offset: ", nexcl+njump,
  55      .      "  lost: ", real(nmax-nmaxp)/real(nmax)*100, " %"
  56          write(iunit,*) "      jump: ", ejump*100, " %"
  57          write(iunit,*) "      slip: ", eslip*100, " %"
  58          write(iunit,*) "  weighted: ", etot*100, " %"
  59          write(iunit,'()')
  60       endif
  61       if(etot.lt.1e-5) stop
  62       nmaxp=nmaxp-1
  63       if(nmaxp.gt.2) goto 1
  64       end
  65
  66       subroutine usage()
  67 c usage message
  68
  69       call whatineed(
  70      .   "[-j# -o outfile -l# -x# -c# -V# -h] file")
  71       call popt("j","weight given to jump relative to slip (0.5)")
  72       call popt("l","number of values to be read (all)")
  73       call popt("x","number of values to be skipped (0)")
  74       call popt("m","number of columns to be read (1)")
  75       call popt("c","columns to be read (1)")
  76       call pout("file_end")
  77       call pall()
  78       stop
  79       end
  80
  81       subroutine jump(nmax,nmaxp,nx,x,mcmax,wjump,etot,ejump,eslip,
  82      .   njump)
  83 c loop through time ofsets to minimize jump effect
  84       dimension x(nx,*)
  85
  86       etot=mcmax
  87       do 10 nj=0,nmax-nmaxp
  88          xj=0
  89          sj=0
  90          do 20 m=1,mcmax
  91             xj=xj+xjump(nmaxp,x(1+nj,m))
  92  20         sj=sj+sjump(nmaxp,x(1+nj,m))
  93          if(wjump*xj+(1-wjump)*sj.ge.etot) goto 10
  94          etot=wjump*xj+(1-wjump)*sj
  95          ejump=xj
  96          eslip=sj
  97          njump=nj
  98  10      continue
  99       end
 100
 101       function xjump(nmax,x)
 102 c contribution of end effect to 1st derivative
 103       dimension x(*)
 104
 105       call rms(nmax,x,sc,sd)
 106       xjump=0
 107       if(sd.eq.0.) return
 108       xjump=(x(1)-x(nmax))**2/(nmax*sd**2)
 109       end
 110
 111       function sjump(nmax,x)
 112 c contribution of end effect to 2nd derivative
 113       dimension x(*)
 114
 115       call rms(nmax,x,sc,sd)
 116       sjump=0
 117       if(sd.eq.0.) return
 118       sjump=((x(nmax)-x(nmax-1))-(x(2)-x(1)))**2 / (nmax*sd**2)
 119       end
 120