← Back to branch summary

~ubuntu-branches/debian/sid/octave-tisean/sid

« back to all changes in this revision

Viewing changes to src/source_f/ts_upo.f

  • Committer: Package Import Robot
  • Author(s): Rafael Laboissiere
  • Date: 2017-08-14 12:53:47 UTC
  • Revision ID: package-import@ubuntu.com-20170814125347-ju5owr4dggr53a2n
Tags: upstream-0.2.3
Import upstream version 0.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/source_f/ts_upo.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   locate unstable periodic points
 
23
c   author T. Schreiber (1998)
 
24
c===========================================================================
 
25
c   modified Piotr Held (2015)
 
26
c===========================================================================
 
27
c==== NOTE MUST BE COMPILED WITH flag -freal-4-real-8 ======
 
28
      subroutine ts_upo(m_in, eps_in, frac,teq,tdis,h,
 
29
     $                  tacc,iper,icen,
 
30
     $                  lines_read, in_out1, olens, orbit_data, sizedat,
 
31
     $                  acc, stability);
 
32
      implicit none
 
33
c -- input declarations --
 
34
      integer*4 m_in, iper, icen, lines_read, sizedat
 
35
      real*8 eps_in, frac, teq, tdis, h, tacc
 
36
c     the *_in is used because they need to be converted from double to float
 
37
 
 
38
      real*8 in_out1(lines_read), orbit_data(sizedat)
 
39
      real*8 acc(icen), stability(icen), olens(icen)
 
40
 
 
41
c -- parameters --
 
42
      integer*4 nx, mper
 
43
      parameter(nx=1000000,mper=20)
 
44
 
 
45
c -- variable declarations --
 
46
 
 
47
c     from main body
 
48
      integer*4 n
 
49
      integer istderr
 
50
      real*8 sc,sd
 
51
      external istderr, rms
 
52
c -- functions --
 
53
      integer known, isold
 
54
      integer*4 iperiod
 
55
      real stab
 
56
 
 
57
c     from old findupo() -- might cause problems
 
58
      integer*4 info, nfev, ndum, itry, ior
 
59
      integer*4 i, ipor
 
60
      integer*4 iw(mper)
 
61
      real tol
 
62
      real enorm, r1mach
 
63
      real xp(mper),fvec(mper),xor(mper,nx)
 
64
      real w0(mper,mper),w1(mper),w2(mper),w3(mper)
 
65
      real w4(mper),w5(mper),w6(mper)
 
66
      real*8 err, sor
 
67
      external peri, enorm, r1mach, snls1
 
68
 
 
69
c     new for creating an output
 
70
      integer*4 orbit_no, data_pos
 
71
 
 
72
 
 
73
c -- global variables --
 
74
      integer*4 nmax, m
 
75
      real*8 eps
 
76
      real*8 x(nx)
 
77
      common /period/ x, nmax, m, eps
 
78
 
 
79
c -- Assign input variables to program variables --
 
80
      nmax=lines_read
 
81
      m=m_in
 
82
      eps=eps_in
 
83
 
 
84
      do 1 i=1,lines_read
 
85
         x(i)=in_out1(i)
 
86
 1       continue
 
87
 
 
88
c -- Prepare the data --
 
89
      call rms(nmax,x,sc,sd)
 
90
      if(frac.gt.0) eps=sd*frac
 
91
      if(teq.lt.0.) teq=eps
 
92
      if(tdis.lt.0.) tdis=eps
 
93
      if(tacc.lt.0.) tacc=eps
 
94
      if(h.lt.0.) h=eps
 
95
 
 
96
c -- Prepare iterators for output --
 
97
      orbit_no = 0;
 
98
      data_pos = 0;
 
99
 
 
100
c -- Exectute main program --
 
101
c     old call findupo(iper,icen,teq,tdis,tacc,h,iunit,iverb)
 
102
      if(iper.gt.mper) then
 
103
         call xstopx ("findupo: make mper larger.")
 
104
      endif
 
105
      tol=sqrt(r1mach(4))
 
106
      itry=0
 
107
      ior=0
 
108
      do 10 n=iper,nmax
 
109
c         if(iv_10(iverb).eq.1) then
 
110
c            if(mod(n,10).eq.0) write(istderr(),'(i7)') n
 
111
c         else if(iv_100(iverb).eq.1) then
 
112
c            if(mod(n,100).eq.0) write(istderr(),'(i7)') n
 
113
c         else if(iv_1000(iverb).eq.1) then
 
114
c            if(mod(n,1000).eq.0) write(istderr(),'(i7)') n
 
115
c         endif
 
116
         if(known(n,iper,teq).eq.1) goto 10
 
117
         itry=itry+1
 
118
         if(itry.gt.icen) goto 999
 
119
         do 20 i=1,iper
 
120
 20         xp(i)=x(n-iper+i)
 
121
         call snls1(peri,1,iper,iper,xp,fvec,w0,mper,tol,tol,0.,
 
122
     .      20*(iper+1),0.,w1,1,100.,0,info,nfev,ndum,iw,w2,w3,w4,w5,w6)
 
123
         err=enorm(iper,fvec)
 
124
         if(info.eq.-1.or.info.eq.5.or.err.gt.tacc) goto 10   ! unsuccessfull
 
125
         if(isold(iper,xp,ior,xor,tdis).eq.1) goto 10         ! already found
 
126
         ior=ior+1                                            ! a new orbit
 
127
         do 30 i=1,iper
 
128
 30         xor(i,ior)=xp(i)
 
129
         ipor=iperiod(iper,xp,tdis)
 
130
         sor=real(ipor)*stab(iper,xp,h)/real(iper)
 
131
 
 
132
c        old call print(iper,xp,ipor,sor,err,iunit,iverb)
 
133
         orbit_no = orbit_no+1
 
134
         olens(orbit_no+1) = ipor
 
135
         acc(orbit_no+1) = err
 
136
         stability(orbit_no+1) = exp(sor)
 
137
         do 40 i=1,ipor
 
138
            data_pos = data_pos + 1
 
139
            orbit_data(data_pos+1) = xp(i)
 
140
 40      continue
 
141
 10   continue
 
142
 
 
143
c -- write lengths input arrays --
 
144
 999  olens(1) = orbit_no
 
145
      orbit_data(1) = data_pos
 
146
      acc(1) = orbit_no
 
147
      stability(1) = orbit_no
 
148
 
 
149
 
 
150
      end
 
151
 
 
152
      integer function known(n,iper,tol)
 
153
c return 1 if equivalent starting point has been tried
 
154
      parameter(nx=1000000)
 
155
      real*8 x(nx)
 
156
      common /period/ x, nmax, m, eps
 
157
 
 
158
      known=1
 
159
      do 10 nn=iper,n-1
 
160
         dis=0
 
161
         do 20 i=1,iper
 
162
 20         dis=dis+(x(n-iper+i)-x(nn-iper+i))**2
 
163
 10      if(sqrt(dis).lt.tol) return
 
164
      known=0
 
165
      end
 
166
 
 
167
      integer function isold(iper,xp,ior,xor,toler)
 
168
c determine if orbit is in data base
 
169
      parameter(mper=20)
 
170
      dimension xp(iper), xor(mper,*)
 
171
 
 
172
      isold=1
 
173
      do 10 ip=1,iper
 
174
         do 20 io=1,ior
 
175
            dor=0
 
176
            do 30 i=1,iper
 
177
 30            dor=dor+(xp(i)-xor(i,io))**2
 
178
 20            if(sqrt(dor).le.toler) return
 
179
 10      call oshift(iper,xp)
 
180
      isold=0
 
181
      end
 
182
  
 
183
      subroutine oshift(iper,xp)
 
184
c leftshift orbit circularly by one position
 
185
      dimension xp(*)
 
186
 
 
187
      h=xp(1)
 
188
      do 10 i=1,iper-1
 
189
 10      xp(i)=xp(i+1)
 
190
      xp(iper)=h
 
191
      end
 
192
 
 
193
      integer*4 function iperiod(iper,xp,tol)
 
194
c determine shortest subperiod
 
195
      dimension xp(*)
 
196
 
 
197
      do 10 iperiod=1,iper
 
198
         dis=0
 
199
         do 20 i=1,iper
 
200
            il=i-iperiod
 
201
            if(il.le.0) il=il+iper
 
202
 20         dis=dis+(xp(i)-xp(il))**2
 
203
 10      if(sqrt(dis).le.tol) return
 
204
      end
 
205
 
 
206
      subroutine peri(iflag,mf,iper,xp,fvec)
 
207
c built discrepancy vector (as called by snls1)
 
208
      dimension xp(iper),fvec(mf)
 
209
 
 
210
      do 10 ip=1,iper
 
211
         fvec(ip)=xp(1)-fc(iper,xp,iflag)
 
212
 10      call oshift(iper,xp)
 
213
      end
 
214
 
 
215
      function fc(iper,xp,iflag)
 
216
c predict (cyclic) point 1, using iper,iper-1...
 
217
      parameter(nx=1000000)
 
218
      dimension  xp(*)
 
219
      real*8 x(nx)
 
220
      common /period/ x, nmax, m, eps
 
221
      data cut/20/
 
222
 
 
223
      eps2=1./(2*eps*eps)
 
224
      ft=0
 
225
      sw=0
 
226
      fc=0
 
227
      do 10 n=m+1,nmax
 
228
         dis=0
 
229
         do 20 i=1,m
 
230
 20         dis=dis+(x(n-i)-xp(mod(m*iper-i,iper)+1))**2
 
231
         ddis=dis*eps2
 
232
         w=0
 
233
         if(ddis.lt.cut) then
 
234
            w=exp(-ddis)
 
235
         endif
 
236
         ft=ft+w*x(n)
 
237
 10      sw=sw+w
 
238
      iflag=-1
 
239
      if(abs(sw).le.1e-312) return   ! fc undefined, stop minimising
 
240
      fc=ft/sw
 
241
      iflag=1
 
242
      end
 
243
 
 
244
      real function stab(ilen,xp,h)
 
245
c compute cycle stability by iteration of a tiny perturbation
 
246
      parameter(nx=1000000,mper=20,maxit=1000)
 
247
      dimension xp(*), xcop(mper)
 
248
      real*8 x(nx)
 
249
      common /period/ x, nmax, m, eps
 
250
 
 
251
      if(mper.lt.ilen) then
 
252
         call xstopx ("stability: make mper larger.")
 
253
      endif
 
254
      iflag=1
 
255
      stab=0
 
256
      do 10 i=2,m
 
257
 10      xcop(i)=xp(mod(i-1,ilen)+1)
 
258
      xcop(1)=xp(1)+h
 
259
      do 20 it=1,maxit
 
260
         do 30 itt=1,ilen
 
261
            xx=fc(m,xcop,iflag)
 
262
            if(iflag.eq.-1) goto 1
 
263
            call oshift(m,xcop)
 
264
 30         xcop(m)=xx
 
265
         dis=0
 
266
         do 40 i=1,m
 
267
 40         dis=dis+(xcop(i)-xp(mod(i-1,ilen)+1))**2
 
268
         dis=sqrt(dis)
 
269
         stab=stab+log(dis/h)
 
270
         do 20 i=1,m
 
271
 20         xcop(i)=xp(mod(i-1,ilen)+1)*(1-h/dis) + xcop(i)*h/dis
 
272
 1    stab=stab/max(it-1,1)
 
273
      end