c c Simulation of non-linear pendulum c - Euler or Runge-Kutta (2nd order usual one) c Fortran version written by H. Nakanishi c Output into files; Also prepare to look at Lyapunov exponent. c program pend5 c c No arrays needed c implicit real*8 (a-h,o-z) character ans,yes character*12 of yes='y' c c Use subroutines to do the work c 20 call initialize(thi,th2i,omi,ext,dt,q,fd,dr,method,of,sep0,np c ,b,c) call calculate(thi,th2i,omi,ext,dt,q,fd,dr,n,method,np,of,b,c) print *,'another try?' read(5,10) ans 10 format(a1) if(ans.eq.yes) go to 20 stop end c subroutine initialize(thi,th2i,omi,ext,dt,q,fd,dr,m,of,sep0,np c ,b,c) c c Initialize variables c implicit real*8 (a-h,o-z) character*12 of print *,'Euler(1), Euler-Cromer(2) or Runge-Kutta2nd (3)? -> ' read(5,*) m if(m.ne.1.and.m.ne.2.and.m.ne.3) then print *,'must select 1, 2, or 3 ..' stop endif print *,'theta0, angular velocity, pendulum length, sep0 ?' read(5,*) thi,omi,ext,sep0 th2i=thi+sep0 print *,'time step, damping const, force amp, ang. frequency?' read(5,*) dt,q,fd,dr print *,'how many (driving) periods?' read(5,*) np print *,'Poincare section? (Yes=>1,No=>0)' read(5,*) b if (b.ne.0) then print *,'initial phase?' read(5,*) c else c=0.d0 endif print *,'output file name?' read(5,14) of 14 format(a12) return end c subroutine calculate(thi,th2i,omi,ext,dt,q,fd,dr,n,m,np,of,b,c) c c Open and rewind the output file c implicit real*8 (a-h,o-z) character*12 of open(1,file=of) rewind(1) c c Now use the Euler method or the Runge-Kutta (2nd order) c g=9.8 period=6.2831853/dr oma=omi omb=omi tha=thi th2a=th2i t=0 e=0.5d0*ext**2*omi**2+g*ext*(1.d0-dcos(thi)) write(1,15) t,tha,th2a,oma,e 15 format(1x,1p,6(e12.5,1x)) 10 continue if(m.eq.1) then call dv(oma,tha,t,g,ext,q,fd,dr,dom,dth) om = oma+dt*dom th = tha+dt*dth call dv(omb,th2a,t,g,ext,q,fd,dr,dom,dth) om2 = omb+dt*dom th2 = th2a+dt*dth elseif(m.eq.2) then call dv(oma,tha,t,g,ext,q,fd,dr,dom,dth) om = oma+dt*dom th = tha+dt*om call dv(omb,th2a,t,g,ext,q,fd,dr,dom,dth) om2 = omb+dt*dom th2 = th2a+dt*om2 else call dv(oma,tha,t,g,ext,q,fd,dr,dom,dth) om1 = oma+0.5d0*dt*dom th1 = tha+0.5d0*dt*dth t1 = t+0.5d0*dt call dv(om1,th1,t1,g,ext,q,fd,dr,dom2,dth2) om = oma+dt*dom2 th = tha+dt*dth2 call dv(omb,th2a,t,g,ext,q,fd,dr,dom,dth) omb1 = omb+0.5d0*dt*dom th21 = th2a+0.5d0*dt*dth t1 = t+0.5d0*dt call dv(omb1,th21,t1,g,ext,q,fd,dr,dom2,dth2) om2 = omb+dt*dom2 th2 = th2a+dt*dth2 endif x1=th-int(th/6.283185307)*6.283185307 if(x1.gt.3.141592654) x1=x1-6.283185307 if(x1.lt.-3.141592654) x1=x1+6.283185307 x2=x1-th2 x2=x2-int(x2/6.283185307)*6.283185307 if(x2.gt.3.141592654) x2=x2-6.283185307 if(x2.lt.-3.141592654) x2=x2+6.283185307 x2=abs(x2) t = t+dt if(b.eq.0) then e=0.5d0*ext**2*om**2+g*ext*(1.d0-dcos(x1)) write(1,15) t,x1,x2,om,e else x3=dr*b*t-c-int((dr*b*t-c)/6.283185307)*6.283185307 if(x3.lt.0) x3=x3+6.283185307 ddt=dr*b*dt*0.5d0 if(x3.lt.ddt .or. x3.gt.(6.283185307-ddt)) then e=0.5d0*ext**2*om**2+g*ext*(1.d0-dcos(x1)) write(1,15) t,x1,x2,om,e endif endif oma = om omb = om2 tha = th th2a = th2 if(t.gt.np*period) then n=i return endif go to 10 close(1) return end c subroutine dv(om0,th0,t0,g,ext,q,fd,dr,dom,dth) implicit real*8 (a-h,o-z) dth=om0 dom=-g/ext*dsin(th0)-q*om0+fd*dsin(dr*t0) return end