--- /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 lorenz.f
+c integrates the Lorenz system with Runge Kutta fourth order
+c author: H. Kantz (2007) based on earlier versions
+c with optional noise
+c===========================================================================
+c
+
+ real*8 x(3),u(3,3),sliap(3),bb,ss,rr,r1,r2,dh,s
+ character*72 fout
+ data iverb/1/
+
+ iverb=ican('V',iverb)
+ call whatido("integration of the Lorenz system",iverb)
+ irun=imust('l')
+ itrans=ican('x',100)
+ rr=fcan('R',28.0)
+ ss=fcan('S',10.0)
+ bb=fcan('B',2.666666667)
+ isamp=ican('f',100)
+ sn=fcan('r',0.)
+c ilyap=lopt('L',1)
+
+ isout=igetout(fout,iverb)
+
+ if(isout.eq.1) fout="lorenz.dat"
+ call outfile(fout,iunit,iverb)
+
+cc intermittency parameters
+c ss=10.d0
+c rr=166.11d0
+c bb=8.d0/3.d0
+
+ iseed1=6456423
+ iseed2=7243431
+
+ xav=0.
+ xsq=0.
+ rsq=0.
+
+c step width of Runge Kutta integration dh:
+ dh=.0005d0
+c time intervals between re-orthogonalization of tangent space
+c vectors: 0.01 units of time.
+ ireno=.01d0/dh
+c length of transient in iteration steps:
+ itrans=real(itrans)/dh
+ totaltime=real(irun)/real(isamp)
+ istep=1.d0/dh/isamp
+
+ if (iverb.eq.1)
+ . write(istderr(),*)'Lorenz trajectory covering',totaltime,
+ . ' time units'
+
+c x(1)=sqrt(s*(r+1.d0))+2.
+c x(2)=x(1)-1.d0
+c x(3)=r
+
+ x(1)=5.
+ x(2)=-10.
+ x(3)=3.
+
+ do 1 i=1,3
+ sliap(i)=0.d0
+ do j=1,3
+ u(i,j)=0.d0
+ enddo
+ u(i,i)=1.d0
+1 continue
+
+ do 10 i2=1,itrans
+
+ call RUKU(3,x,u,dh,bb,ss,rr)
+
+ if (mod(i2,ireno).eq.0) then
+ call norm(u,1,s)
+ do i=2,3
+ do j=1,i-1
+ call proj(u,i,j)
+ enddo
+ call NORM(u,i,s)
+ enddo
+ endif
+
+10 continue
+
+ write(iunit,101)x(1),x(2),x(3)
+
+ 100 continue
+ do 20 i2=1,irun*istep
+c add noise
+ if (sn.gt.0.0) then
+ call gauss(r1,r2,iseed1,iseed2)
+ x(1)=x(1)+r1*sn
+ x(2)=x(2)+r2*sn
+ call gauss(r1,r2,iseed1,iseed2)
+ x(3)=x(3)+r1*sn
+ xav=xav+x(1)
+ xsq=xsq+x(1)**2
+ rsq=rsq+r1*r1
+ endif
+ call RUKU(3,x,u,dh,bb,ss,rr)
+ if (mod(i2,istep).eq.0) write(iunit,101)x(1),x(2),x(3)
+ if (mod(i2,ireno).eq.0) then
+c Gram Schmidt Orthonormierung
+ call norm(u,1,s)
+ sliap(1)=sliap(1)+log(s)
+ do i=2,3
+ do j=1,i-1
+ call proj(u,i,j)
+ enddo
+ call NORM(u,i,s)
+ sliap(i)=sliap(i)+log(s)
+ enddo
+ endif
+
+ 20 continue
+
+ if (sn.gt.0.0) then
+ xav=xav/irun/istep
+ xsq=xsq/irun/istep
+ rsq=rsq/irun/istep
+ rlevel=sqrt(rsq)/sqrt(xsq-xav*xav)*100.
+ if (iverb.eq.1)
+ . print*,'noise level in percent of x-coordinate',rlevel
+ endif
+ if (iverb.eq.1) then
+ write(istderr(),*)
+ write(istderr(),*)'Lyapunov exponents [1/unit time]'
+ do i=1,3
+ write(istderr(),*)real(sliap(i)/totaltime)
+ enddo
+ endif
+
+ 101 format(2x,3f10.3)
+
+ stop
+ end
+
+ subroutine FORCE(x,ff,dh,bb,ss,rr)
+ real*8 x(3),ff(3),dh,bb,ss,rr
+
+ ff(1)=dh*ss*(x(2)-x(1))
+ ff(2)=dh*(x(1)*(-x(3)+rr)-x(2))
+ ff(3)=dh*(x(1)*x(2)-bb*x(3))
+
+ return
+ end
+
+ subroutine LFORCE(x,u,fl,dh,bb,ss,rr)
+ real*8 x(3),u(3,3),dh,fl(3,3),bb,ss,rr
+
+ do j=1,3
+ fl(j,1)=dh*ss*(u(j,2)-u(j,1))
+ fl(j,2)=dh*(u(j,1)*(rr-x(3))-x(1)*u(j,3)-u(j,2))
+ fl(j,3)=dh*(u(j,1)*x(2)+x(1)*u(j,2)-bb*u(j,3))
+ enddo
+ return
+ end
+
+ subroutine RUKU(n,x,u,dh,bb,ss,rr)
+c 4th-order Runge Kutta
+c initial point x
+c final point y
+c stepsize dh
+c add subroutine force
+
+ implicit real*8 (a-h,o-z)
+ real*8 x(3),ff1(3),ff2(3),ff3(3),ff4(3),dummy(3)
+ real*8 u(3,3),fl1(3,3),fl2(3,3),fl3(3,3),fl4(3,3)
+ real*8 dl(3,3)
+
+ call force(x,ff1,dh,bb,ss,rr)
+ call LFORCE(x,u,fl1,dh,bb,ss,rr)
+
+ do i=1,n
+ dummy(i)=ff1(i)*.5d0+x(i)
+ do j=1,3
+ dl(i,j)=fl1(i,j)*.5d0+u(i,j)
+ enddo
+ enddo
+
+ call force(dummy,ff2,dh,bb,ss,rr)
+ call LFORCE(dummy,dl,fl2,dh,bb,ss,rr)
+
+ do i=1,n
+ dummy(i)=ff2(i)*.5d0+x(i)
+ do j=1,3
+ dl(i,j)=fl2(i,j)*.5d0+u(i,j)
+ enddo
+ enddo
+
+ call force(dummy,ff3,dh,bb,ss,rr)
+ call LFORCE(dummy,dl,fl3,dh,bb,ss,rr)
+
+ do i=1,n
+ dummy(i)=ff3(i)+x(i)
+ do j=1,3
+ dl(i,j)=fl3(i,j)+u(i,j)
+ enddo
+ enddo
+
+ call force(dummy,ff4,dh,bb,ss,rr)
+ call LFORCE(dummy,dl,fl4,dh,bb,ss,rr)
+
+ do i=1,n
+ x(i)=x(i)+ff1(i)/6.d0+ff2(i)/3.d0+ff3(i)/3.d0+ff4(i)/6.d0
+ do j=1,3
+ u(i,j)=u(i,j)+fl1(i,j)/6.d0+fl2(i,j)/3.d0+fl3(i,j)/3.d0
+ + +fl4(i,j)/6.d0
+ enddo
+ enddo
+
+ return
+ end
+
+ subroutine NORM(u,i,s)
+ real*8 u(3,3),s
+
+ s=0.d0
+ do 10 j=1,3
+10 s=s+u(i,j)**2
+ s=sqrt(s)
+ si=1.d0/s
+ do 20 j=1,3
+20 u(i,j)=u(i,j)*si
+ return
+ end
+
+ subroutine PROJ(u,i,j)
+ real*8 u(3,3),s
+ s=0.d0
+ do 10 k=1,3
+10 s=s+u(i,k)*u(j,k)
+ do 20 k=1,3
+20 u(i,k)=u(i,k)-s*u(j,k)
+ return
+ end
+
+c>-------------------------------------------------------
+ subroutine gauss(r1,r2,iseed1,iseed2)
+
+ real*8 r1,r2,p,phi,r
+ pii=8.d0*atan(1.d0)
+
+ call RANDOM1(p,iseed1)
+ call RANDOM1(phi,iseed2)
+
+ phi=phi*pii
+ r=sqrt(-log(1.d0-p)*2.d0)
+
+ r1=r*sin(phi)
+ r2=r*cos(phi)
+ return
+ end
+c>-------------------------------------------------------
+ subroutine RANDOM1(r,iseed)
+c
+c random number generator of Park & Miller
+c random numbers in [0,1] !!!
+ real*8 r
+ integer*8 ia,im,ix
+ ia=7**5
+ im=2147483647
+ ix=iseed
+ ix=mod(ia*ix,im)
+ r=dfloat(ix)/dfloat(im)
+ iseed=ix
+ return
+ end
+c>------------------------------------------------------------------
+ subroutine usage()
+c usage message
+
+ call whatineed(
+ . "-l# [-f# -r# -R# -S# -B# -o outfile -x# -V# -h]")
+ call popt("l","length of trajectory x,y,z")
+ call popt("f","sample points per unit time [100]")
+ call popt("r","absolute noise amplitute [0]")
+ call popt("R","parameter r [28]")
+ call popt("S","parameter sigma [10]")
+ call popt("B","parameter b [8/3]")
+ call popt("x","transient discarded [100 units of time]")
+c call popt("L","if present: compute Lyapunov exponents")
+ call pout("lorenz.dat")
+ call pall()
+ stop
+ end
+
+
+
git://source.jalview.org / jabaws.git / blobdiff