c
c Simulation of Poincare and bifarcation 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 bi
c
        character ans,yes
        character*12 of
	real*8 th0,om0,ext,dt,q,fmin,fmax,df,dr,phi,x
        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
	character*12 of
	real*8 th0,om0,ext,dt,q,fmin,fmax,df,dr,phi,x
	print *,'theta0, angular velocity, pendulum length ?'
	read(5,*) th0,om0,ext
	print *,'time step, damping const, driving angular frequency?'
	read(5,*) dt,q,dr
	print *,'fdmin, fdmax, df ?'
	read(5,*) fmin,fmax,df
	print *,'sec 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,g,period
	real*8 t,th,om,t1,th1,om1,dom,dth,dom2,dth2,t3,x1,e,fd
	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 om0,th0,t0,g,ext,q,fd,dr,dom,dth
	dth=om0
	dom=-g/ext*dsin(th0)-q*om0+fd*dsin(dr*t0)
	return
	end
