--- /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