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 Line 74 originally had incorrrect index; corrected to: "t1 = t(i-1)+0.5*dt"
c
	program pendulum
c
c       Declare the arrays we will need
c
	dimension th(5003),om(5003),t(5003)
c
c       Use subroutines to do the work
c
	call initialize(th,om,t,ext,dt,q,fd,dr,lsym,nsym,method)
	call calculate(th,om,t,ext,dt,q,fd,dr,n,lsym,nsym,method)
	call display(th,om,t,ext,dt,q,fd,dr,n,lsym,nsym,method)
	stop
	end
c
	subroutine initialize(th,om,t,ext,dt,q,fd,dr,lsym,nsym,method)
c
c       Initialize variables
c
        dimension th(1),om(1),t(1)
	character ans,yes
	yes='y'
	print *,'Euler(1), Euler-Cromer(2) or Runge-Kutta(3)? -> '
	read(5,*) method
	if(method.ne.1.and.method.ne.2.and.method.ne.3) then
		print *,'must select 1, 2, or 3 ..'
		stop
	endif
	print *,'initial angle, angular velocity, pendulum length?'
	read(5,*) th(1),om(1),ext
	t(1)=0
	print *,'time step, damping const, force amp, frequency?'
	read(5,*) dt,q,fd,dr
	print *,'set line, symbol?'
	read(5,14) ans
14	format(a1)
	if(ans.eq.yes) then
		print *,'lsym, nsym ? -> '
		read(5,*) lsym,nsym
	else
		lsym=-1
		nsym=1
	endif
	return
	end
c
	subroutine calculate(th,om,t,ext,dt,q,fd,dr,n,lsym,nsym,method)
c
c       Now use the Euler method or the Runge-Kutta (2nd order)
c
	dimension th(1),om(1),t(1)
	g=9.8
	period=6.2831853/sqrt(g/ext)
	nmax=5000
	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
	    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)
	    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-1)+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
	    endif
	    if(t(i).ge.10*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 display(th,om,t,ext,dt,q,fd,dr,n,lsym,nsym,method)
c
c       First set up title and label axes for graph. Plotting is a lot
c       of work in fortran.
c       This version displays output as well as writes to a file "graph.out".
c
	dimension th(1),om(1),t(1)
	call usrmon(.true.,.false.,-0.5,9.,-0.5,12.)
	call pltlun(19,.true.,.false.)
	call pltlfn('graph.out')
	call plots
	call plot(0.,3.5,-3)
	if(method.eq.1) then
	call text(0.2,5.8,0.2,'Pendulum: Euler',0.,15,0)
	elseif(method.eq.2) then
	call text(0.2,5.8,0.2,'Pendulum: Euler-Cromer',0.,22,0)
	else
	call text(0.2,5.8,0.2,'Pendulum: Runge-Kutta2',0.,22,0)
	endif
	call text(0.2,5.5,0.2,'Angle vs Time',0.,13,0)
	call scalex(t,5.,n,1)
	call axctl(0.15,0.04,0.15,0.2,0.2,-1)
	call axisx(0.,0.,'Time (s)',-8,5.,0.,
     c  t(n+1),t(n+2),t(n+3),4)
	call axisx(0.,5.,' ',1,5.,0.,t(n+1),t(n+2),t(n+3),20)
	call scalex(th,5.,n,1)
	call axctl(0.15,0.04,0.15,0.2,0.2,2)
	call axisx(0.,0.,'Angle (radians)',15,5.,90.,
     c	th(n+1),th(n+2),th(n+3),-5)
	call axisx(5.,0.,' ',-1,5.,90.,
     c  th(n+1),th(n+2),th(n+3),-21)
	call line(t,th,n,1,lsym,nsym)
	call number(0.5,3.7,0.1,q,0.,'''q = '',f5.2')
	call number(0.5,3.9,0.1,fd,0.,'''fd= '',f5.2')
	call number(0.5,4.1,0.1,dr,0.,'''dr= '',f8.5')
	call number(0.5,4.3,0.1,ext,0.,'''ext = '',f5.2')
	call number(0.5,4.5,0.1,dt,0.,'''dt = '',f5.2')
	call number(0.5,4.7,0.1,ext,0.,'''ext = '',f5.2')
	call plot(0.,0.,999)
	return
	end