c
c Simulation of non-linear pendulum
c        - Euler or Runge-Kutta (2nd order usual one)
c Fortran version written by H. Nakanishi, need to be compiled with "-lpepl".
c Output into files; Also prepare to look at Lyapunov exponent.
c
	program pendulum3
c
c       Declare the arrays we will need
c
	dimension th(65536),om(65536),t(65536),th2(65536)
        character ans,yes
        character*12 of
        yes='y'
c
c       Use subroutines to do the work
c
20	call initialize(th,th2,om,t,ext,dt,q,fd,dr,method,of,sep0
     c  ,np,ishift)
	call calculate(th,th2,om,t,ext,dt,q,fd,dr,n,method,np)
	call record(th,th2,om,t,n,of,ishift)
        print *,'another try?'
        read(5,10) ans
10      format(a1)
        if(ans.eq.yes) go to 20
	stop
	end
c
	subroutine initialize(th,th2,om,t,ext,dt,q,fd,dr,m,of,sep0
     c  ,np,ishift)
c
c       Initialize variables
c
        dimension th(1),th2(1),om(1),t(1)
	character*12 of
	character*1 ans,yes
	yes="y"
	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,*) th(1),om(1),ext,sep0
	t(1)=0
	th2(1)=th(1)+sep0
	print *,'time step, damping const, force amp, frequency?'
	read(5,*) dt,q,fd,dr
	print *,'how many (natural) periods?'
	read(5,*) np
	print *,'shift theta to (-pi,pi]?'
        read(5,15) ans
15      format(a1)
        if(ans.eq.yes) then
		ishift=1
	else
		ishift=0
	endif
        print *,'output file name?'
        read(5,14) of
14      format(a12)
	return
	end
c
	subroutine calculate(th,th2,om,t,ext,dt,q,fd,dr,n,method,np)
c
c       Now use the Euler method or the Runge-Kutta (2nd order)
c
	dimension th(1),th2(1),om(1),t(1)
	g=9.8
	period=6.2831853/sqrt(g/ext)
	nmax=min(65536,int(np*period/dt))
	oma=om(1)
	do 10 i = 2,nmax
		t(i) = t(i-1)+dt
	    if(method.eq.1) then
		call dv(om(i-1),th(i-1),t(i-1),g,ext,q,fd,dr,dom,dth)
         	om(i) = om(i-1)+dt*dom
         	th(i) = th(i-1)+dt*dth

		call dv(oma,th2(i-1),t(i-1),g,ext,q,fd,dr,dom,dth)
         	oma = oma+dt*dom
         	th2(i) = th2(i-1)+dt*dth
	    elseif(method.eq.2) then
		call dv(om(i-1),th(i-1),t(i-1),g,ext,q,fd,dr,dom,dth)
         	om(i) = om(i-1)+dt*dom
         	th(i) = th(i-1)+dt*om(i)

		call dv(oma,th2(i-1),t(i-1),g,ext,q,fd,dr,dom,dth)
         	oma = oma+dt*dom
         	th2(i) = th2(i-1)+dt*oma
	    else
		call dv(om(i-1),th(i-1),t(i-1),g,ext,q,fd,dr,dom,dth)
         	om1 = om(i-1)+0.5*dt*dom
         	th1 = th(i-1)+0.5*dt*dth
		t1 = t(i)+0.5*dt
		call dv(om1,th1,t1,g,ext,q,fd,dr,dom2,dth2)
		om(i)=om(i-1)+dt*dom2
		th(i)=th(i-1)+dt*dth2

		call dv(oma,th2(i-1),t(i-1),g,ext,q,fd,dr,dom,dth)
         	oma1 = oma+0.5*dt*dom
         	th21 = th2(i-1)+0.5*dt*dth
		t1 = t(i)+0.5*dt
		call dv(oma1,th21,t1,g,ext,q,fd,dr,dom2,dth2)
		oma=oma+dt*dom2
		th2(i)=th2(i-1)+dt*dth2
	    endif
	    if(t(i).ge.np*period) then
		n=i
		return
	    endif
10	continue
  	n=nmax
	return
	end
c
	subroutine dv(om0,th0,t0,g,ext,q,fd,dr,dom,dth)
	dth=om0
	dom=-g/ext*sin(th0)-q*om0+fd*sin(dr*t0)
	return
	end
c
	subroutine record(th,th2,om,t,n,ofile,ishift)
c
        character*12 ofile
        dimension th(1),th2(1),om(1),t(1)
        open(1,file=ofile)
        do 10 i=1,n
	x1=th(i)-int(th(i)/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(i)
	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)
        if(ishift.eq.1) then
          write(1,15) t(i),x1,x2,om(i)
        else
          write(1,15) t(i),th(i),x2,om(i)
        endif
10      continue
15      format(1x,1p,4(e12.5,",",1x))
        close(1)
	return
	end

