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

