c c Simulation of Poincare and bifurcation in non-linear pendulum c - Runge-Kutta (2nd order usual one) c Fortran version written by H. Nakanishi c Output into files; Poincare, and Bifarcation c program bifurc c real*8 th0,om0,ext,dt,q,fmin,fmax,df,dr,phi,x character ans,yes character*12 of yes='y' c c Use subroutines to do the work c 20 call init(th0,om0,ext,dt,q,fmin,fmax,df,dr,of,i1,i2,phi,x) call calc(th0,om0,ext,dt,q,fmin,fmax,df,dr,of,i1,i2,phi,x) print *,'another try?' read(5,10) ans 10 format(a1) if(ans.eq.yes) go to 20 stop end c subroutine init(th0,om0,ext,dt,q,fmin,fmax,df,dr,of,i1,i2,phi,x) c real*8 th0,om0,ext,dt,q,fmin,fmax,df,dr,phi,x character*12 of print *,'theta0, angular velocity, pendulum length ?' read(5,*) th0,om0,ext print *,'time step, damping const, driving frequency?' read(5,*) dt,q,dr print *,'fdmin, fdmax, df ?' read(5,*) fmin,fmax,df print *,'omega, begin and end periods, phase for section ?' read(5,*) x,i1,i2,phi print *,'output file name ?' read(5,14) of 14 format(a12) return end c subroutine calc(th0,om0,ext,dt,q,fmin,fmax,df,dr,of,i1,i2,phi,x) c real*8 th0,om0,ext,dt,q,fmin,fmax,df,dr,phi,x,x1,e real*8 g,period,t,th,om,th1,dth,dth2,om1,dom,dom2,t1,t3,ddt character*12 of open(1,file=of) rewind 1 if(df.ne.0) then mmax=int((fmax-fmin)/df)+1 else mmax=1 endif g=9.8d0 period=6.2831853/x 15 format(1x,1p,5(e12.5,1x)) do 20 j=1,mmax fd=fmin+(j-1)*df t=0.d0 th=th0 om=om0 10 continue call dv(om,th,t,g,ext,q,fd,dr,dom,dth) om1=om+0.5*dt*dom th1=th+0.5*dt*dth t1=t+0.5*dt call dv(om1,th1,t1,g,ext,q,fd,dr,dom2,dth2) om=om+dt*dom2 th=th+dt*dth2 t=t+dt if(t.lt.i1*period) go to 10 if(t.gt.i2*period) go to 20 t3=x*t-phi-int((x*t-phi)/6.283185307)*6.283185307 if(t3.lt.0) t3=t3+6.283185307 ddt=0.5*dt*x if(t3.lt.ddt.or.t3.gt.(6.283185307-ddt)) then 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 e=0.5*ext**2*om**2+g*ext*(1.-dcos(x1)) write(1,15) t,x1,om,e,fd endif go to 10 20 continue close(1) return end c subroutine dv(om0,th0,t0,g,ext,q,fd,dr,dom,dth) real*8 th0,om0,ext,dt,q,fmin,fmax,df,dr,phi,x real*8 dth,dom,g,t0 dth=om0 dom=-g/ext*dsin(th0)-q*om0+fd*dsin(dr*t0) return end