+++ /dev/null
-c===========================================================================
-c
-c This file is part of TISEAN
-c
-c Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
-c
-c TISEAN is free software; you can redistribute it and/or modify
-c it under the terms of the GNU General Public License as published by
-c the Free Software Foundation; either version 2 of the License, or
-c (at your option) any later version.
-c
-c TISEAN is distributed in the hope that it will be useful,
-c but WITHOUT ANY WARRANTY; without even the implied warranty of
-c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-c GNU General Public License for more details.
-c
-c You should have received a copy of the GNU General Public License
-c along with TISEAN; if not, write to the Free Software
-c Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-c
-c===========================================================================
-c locate unstable periodic points
-c author T. Schreiber (1998)
-c===========================================================================
- parameter(nx=1000000,mper=20)
- dimension x(nx)
- character*72 file, fout
- common /period/ x, nmax, m, eps
- data frac/0./, iper/1/, teq/-1./, tdis/-1./, tacc/-1./, h/-1./
- data iverb/1/
-
- call whatido("locate unstable periodic points",iverb)
- m=max(imust("m"),1)
- eps=fcan("r",0.)
- frac=fcan("v",frac)
- teq=fcan("w",teq)
- tdis=fcan("W",tdis)
- h=fcan("s",h)
- tacc=fcan("a",tacc)
- iper=ican("p",iper)
- nmaxx=ican("l",nx)
- nexcl=ican("x",0)
- jcol=ican("c",0)
- icen=ican("n",nmaxx)
- isout=igetout(fout,iverb)
- if(eps.eq.0.and.frac.eq.0.) call usage()
-
- do 10 ifi=1,nstrings()
- call nthstring(ifi,file)
- nmax=nmaxx
- call readfile(nmax,x,nexcl,jcol,file,iverb)
- if(file.eq."-") file="stdin"
- if(isout.eq.1) call addsuff(fout,file,"_upo_")
- call rms(nmax,x,sc,sd)
- if(frac.gt.0) eps=sd*frac
- if(teq.lt.0.) teq=eps
- if(tdis.lt.0.) tdis=eps
- if(tacc.lt.0.) tacc=eps
- if(h.lt.0.) h=eps
- if(isout.eq.1)
- . write(fout(index(fout,"_upo_")+5:72),'(i2.2)') iper
- call outfile(fout,iunit,iverb)
- call findupo(iper,icen,teq,tdis,tacc,h,iunit,iverb)
- 10 if(iunit.ne.istdout()) close(iunit)
- end
-
- subroutine usage()
-c usage message
- call whatineed(
- . "-m# [-r# | -v#] [-p# -w# -W# -a# -s# -n#"//
- . " -o outfile -l# -x# -c# -V# -h] file(s)")
- call ptext("either -r or -v must be present")
- call popt("m","embedding dimension")
- call popt("r","absolute kernel bandwidth")
- call popt("v","same as fraction of standard deviation")
- call popt("p","period of orbit (1)")
- call popt("w","minimal separation of trial points (e)")
- call popt("W","minimal separation of distinct orbits (e)")
- call popt("a",
- . "maximal error of orbit to be plotted (all plotted)")
- call popt("s","initial separation for stability (e)")
- call popt("n","number of trials (all points)")
- call popt("l","number of values to be read (all)")
- call popt("x","number of values to be skipped (0)")
- call popt("c","column to be read (1 or file,#)")
- call pout("file_upo_pp")
- call pall()
- call ptext("Verbosity levels (add what you want):")
- call ptext(" 1 = input/output" )
- call ptext(" 2 = print orbits found")
- call ptext(" 4 = status after 1000 points")
- call ptext(" 8 = status after 100 points")
- call ptext(" 16 = status after 10 points")
- stop
- end
-
- subroutine findupo(iper,icen,teq,tdis,tacc,h,iunit,iverb)
- parameter(nx=1000000,mper=20)
- external peri
- dimension x(nx),xp(mper),fvec(mper),xor(mper,nx),
- . iw(mper),w0(mper,mper),w1(mper),w2(mper),w3(mper),
- . w4(mper),w5(mper),w6(mper)
- common /period/ x, nmax, m, eps
-
- if(iper.gt.mper) stop "findupo: make mper larger."
- tol=sqrt(r1mach(4))
- itry=0
- ior=0
- do 10 n=iper,nmax
- if(iv_10(iverb).eq.1) then
- if(mod(n,10).eq.0) write(istderr(),'(i7)') n
- else if(iv_100(iverb).eq.1) then
- if(mod(n,100).eq.0) write(istderr(),'(i7)') n
- else if(iv_1000(iverb).eq.1) then
- if(mod(n,1000).eq.0) write(istderr(),'(i7)') n
- endif
- if(known(n,iper,teq).eq.1) goto 10
- itry=itry+1
- if(itry.gt.icen) return
- do 20 i=1,iper
- 20 xp(i)=x(n-iper+i)
- call snls1(peri,1,iper,iper,xp,fvec,w0,mper,tol,tol,0.,
- . 20*(iper+1),0.,w1,1,100.,0,info,nfev,ndum,iw,w2,w3,w4,w5,w6)
- err=enorm(iper,fvec)
- if(info.eq.-1.or.info.eq.5.or.err.gt.tacc) goto 10 ! unsuccessfull
- if(isold(iper,xp,ior,xor,tdis).eq.1) goto 10 ! already found
- ior=ior+1 ! a new orbit
- do 30 i=1,iper
- 30 xor(i,ior)=xp(i)
- ipor=iperiod(iper,xp,tdis)
- sor=ipor*stab(iper,xp,h)/real(iper)
- call print(iper,xp,ipor,sor,err,iunit,iverb)
- 10 continue
- end
-
- function known(n,iper,tol)
-c return 1 if equivalent starting point has been tried
- parameter(nx=1000000)
- dimension x(nx)
- common /period/ x, nmax, m, eps
-
- known=1
- do 10 nn=iper,n-1
- dis=0
- do 20 i=1,iper
- 20 dis=dis+(x(n-iper+i)-x(nn-iper+i))**2
- 10 if(sqrt(dis).lt.tol) return
- known=0
- end
-
- function isold(iper,xp,ior,xor,toler)
-c determine if orbit is in data base
- parameter(mper=20)
- dimension xp(iper), xor(mper,*)
-
- isold=1
- do 10 ip=1,iper
- do 20 io=1,ior
- dor=0
- do 30 i=1,iper
- 30 dor=dor+(xp(i)-xor(i,io))**2
- 20 if(sqrt(dor).le.toler) return
- 10 call oshift(iper,xp)
- isold=0
- end
-
- subroutine oshift(iper,xp)
-c leftshift orbit circularly by one position
- dimension xp(*)
-
- h=xp(1)
- do 10 i=1,iper-1
- 10 xp(i)=xp(i+1)
- xp(iper)=h
- end
-
- function iperiod(iper,xp,tol)
-c determine shortest subperiod
- dimension xp(*)
-
- do 10 iperiod=1,iper
- dis=0
- do 20 i=1,iper
- il=i-iperiod
- if(il.le.0) il=il+iper
- 20 dis=dis+(xp(i)-xp(il))**2
- 10 if(sqrt(dis).le.tol) return
- end
-
- subroutine peri(iflag,mf,iper,xp,fvec,fjac,ldfjac)
-c built discrepancy vector (as called by snls1)
- dimension xp(*),fvec(*)
-
- do 10 ip=1,iper
- fvec(ip)=xp(1)-fc(iper,xp,iflag)
- 10 call oshift(iper,xp)
- end
-
- function fc(iper,xp,iflag)
-c predict (cyclic) point 1, using iper,iper-1...
- parameter(nx=1000000)
- dimension xp(*), x(nx)
- common /period/ x, nmax, m, eps
- data cut/20/
-
- eps2=1./(2*eps*eps)
- ft=0
- sw=0
- fc=0
- do 10 n=m+1,nmax
- dis=0
- do 20 i=1,m
- 20 dis=dis+(x(n-i)-xp(mod(m*iper-i,iper)+1))**2
- ddis=dis*eps2
- w=0
- if(ddis.lt.cut) w=exp(-ddis)
- ft=ft+w*x(n)
- 10 sw=sw+w
- iflag=-1
- if(sw.eq.0) return ! fc undefined, stop minimising
- fc=ft/sw
- iflag=1
- end
-
- function stab(ilen,xp,h)
-c compute cycle stability by iteration of a tiny perturbation
- parameter(nx=1000000,mper=20,maxit=1000)
- dimension xp(*), x(nx), xcop(mper)
- common /period/ x, nmax, m, eps
-
- if(mper.lt.ilen) stop "stability: make mper larger."
- iflag=1
- stab=0
- do 10 i=2,m
- 10 xcop(i)=xp(mod(i-1,ilen)+1)
- xcop(1)=xp(1)+h
- do 20 it=1,maxit
- do 30 itt=1,ilen
- xx=fc(m,xcop,iflag)
- if(iflag.eq.-1) goto 1
- call oshift(m,xcop)
- 30 xcop(m)=xx
- dis=0
- do 40 i=1,m
- 40 dis=dis+(xcop(i)-xp(mod(i-1,ilen)+1))**2
- dis=sqrt(dis)
- stab=stab+log(dis/h)
- do 20 i=1,m
- 20 xcop(i)=xp(mod(i-1,ilen)+1)*(1-h/dis) + xcop(i)*h/dis
- 1 stab=stab/max(it-1,1)
- end
-
- subroutine print(iper,xp,ipor,sor,err,iunit,iverb)
-c write orbit to iunit and to stdout
- dimension xp(*)
-
- write(iunit,*)
- write(iunit,*) "period / accuracy / stability"
- write(iunit,*) ipor, err, exp(sor)
- do 10 i=1,ipor
- 10 write(iunit,*) i, xp(i)
- if(iv_upo(iverb).eq.0) return
- write(istderr(),*)
- write(istderr(),*) "period / accuracy / stability"
- write(istderr(),*) ipor, err, exp(sor)
- do 20 i=1,ipor
- 20 write(istderr(),*) i, xp(i)
- end
git://source.jalview.org / jabaws.git / blobdiff