*     PROGRAM poincare
*     plot phase diagram and Poincare map for damped, driven pendulum
      IMPLICIT NONE
      REAL*8 A, ang_vel, dA, dt, gamma, t, theta
      INTEGER IC1_, IC2_, IR_, istep, nstep
      CHARACTER*80 flag
      LOGICAL Ltmp1_, TBkeyinput
      DATA IC1_/1/
      DATA IC2_/2/
      CALL GWopen(IR_, 0)
      CALL initial(theta,ang_vel,gamma,A,t,dt,nstep,dA,IC1_,IC2_,flag)
      CALL set_up_windows(A,IC1_,IC2_)
      CALL draw_axes(A,IC1_,IC2_)
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         DO istep = 1, nstep
            CALL RK4(theta,ang_vel,gamma,A,t,dt)
            CALL phase_space(theta,ang_vel,IC1_)
         END DO
         CALL Poincare(theta,ang_vel,IC2_)
         IF(TBkeyinput()) CALL change_amplitude(A,dA,IC1_,IC2_,flag)
         Ltmp1_ = flag .NE. 'stop'
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(theta0,ang_vel0,gamma,A,t0,dt,nstep,dA,IC1_
     +,IC2_,flag)
      IMPLICIT NONE
      REAL*8 T, pi, theta0, ang_vel0, gamma, A, t0, dt, dA
      REAL GWaspect
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER IR_, nstep, IC1_, IC2_
      CHARACTER*80 flag
      WRITE(*,'(A,$)') 'initial angle = '
      READ(*,*) theta0
      WRITE(*,'(A,$)') 'initial angular velocity = '
      READ(*,*) ang_vel0
      t0 = 0                                ! initial time
      WRITE(*,'(A,$)') 'damping constant = '
      READ(*,*) gamma
      WRITE(*,'(A,$)') 'amplitude of external force = '
      READ(*,*) A
*     angular frequency of external force equals two and
*     hence period of external force equals pi
      nstep = 100                           ! # iterations between Poincare plot
      T = pi                                ! period of external force
      dt = T/nstep                          ! time step
      dA = 0.01D0                           ! increase in amplitude of external force
*     left-side of screen for phase space plot
      CALL GWvport(IR_, 0.0, 0.0, 0.49*GWaspect(1), 1.0)
*     open second window for Poincare map
      CALL GWsavevp(IR_, IC1_)
      CALL GWsetogn(IR_, 0, IC1_)
      CALL GWvport(IR_, 0.51*GWaspect(1), 0.0, GWaspect(1), 1.0)
      flag = ''
      CALL GWsavevp(IR_, IC2_)
      CALL GWsetogn(IR_, 0, IC2_)
      END
C
      SUBROUTINE set_up_windows(A,IC1_,IC2_)
      IMPLICIT NONE
      REAL*8 A
      REAL vmax, pi
      PARAMETER(pi = 3.141592)
      INTEGER IR_, IC1_, IC2_
      vmax = REAL(4*A)
      CALL GWselvp(IR_, IC1_)
      CALL GWsetogn(IR_, 0, IC1_)
      CALL GWindow(IR_, -pi, -vmax, pi, vmax)
      CALL GWselvp(IR_, IC2_)
      CALL GWsetogn(IR_, 0, IC2_)
      CALL GWindow(IR_, -pi, -vmax, pi, vmax)
      END
C
      SUBROUTINE draw_axes(A,IC1_,IC2_)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 A
      REAL vmax, pi
      PARAMETER(pi = 3.141592)
      INTEGER IR_, IC1_, IC2_
      CALL GWclear(IR_, -1)
      CALL GWselvp(IR_, IC1_)
      CALL GWsetogn(IR_, 0, IC1_)
      CALL GWerase(IR_, IC1_, -1)
      CALL GWsetpen(IR_, 0, -1, -1, -1)
      vmax = 4*A
      CALL GWline(IR_, 0.0, -vmax, 0.0, vmax)
      CALL GWline(IR_, -pi, 0.0, pi, 0.0)
      WRITE(Stmp1_, '(A,F7.5)') 'A = ',A
      CALL GWputtxt(IR_, -3.0, vmax*22.5/24, Stmp1_)
      CALL GWselvp(IR_, IC2_)
      CALL GWsetogn(IR_, 0, IC2_)
      CALL GWerase(IR_, IC2_, -1)
      CALL GWsetpen(IR_, 0, -1, -1, -1)
      CALL GWline(IR_, -pi, 0.0, pi, 0.0)
      CALL GWline(IR_, 0.0, -vmax, 0.0, vmax)
      END
C
      SUBROUTINE RK4(theta,ang_vel,gamma,A,t,dt)
*     fourth-order Runge Kutta algorithm
      IMPLICIT NONE
      REAL*8 k1v, k1x, k2v, k2x, k3v, k3x, k4v, k4x, pi, theta,
     +       ang_vel, gamma, A, t, dt, force
      PARAMETER(pi = 3.141592653589793D0)
      k1v = force(theta,ang_vel,gamma,A,t)*dt
      k1x = ang_vel*dt
      t = t + 0.5D0*dt
      k2v = force(theta+0.5D0*k1x,ang_vel+0.5D0*k1v,gamma,A,t)*dt
      k2x = (ang_vel + 0.5D0*k1v)*dt
      k3v = force(theta+0.5D0*k2x,ang_vel+0.5D0*k2v,gamma,A,t)*dt
      k3x = (ang_vel + 0.5D0*k2v)*dt
      t = t + 0.5D0*dt
      k4v = force(theta+k3x,ang_vel+k3v,gamma,A,t)*dt
      k4x = (ang_vel + k3v)*dt
      ang_vel = ang_vel + (k1v + 2*k2v + 2*k3v + k4v)/6
      theta = theta + (k1x + 2*k2x + 2*k3x + k4x)/6
      IF(theta .GT. pi) theta = theta - 2*pi
      IF(theta .LT. -pi) theta = theta + 2*pi
      END
C
      FUNCTION force(theta,ang_vel,gamma,A,t)
*     A is amplitude of driving force
      IMPLICIT NONE
      REAL*8 force, theta, ang_vel, gamma, A, t
      force = -gamma*ang_vel - (1 +2*A*cos(2*t))*sin(theta)
      END
C
      SUBROUTINE phase_space(theta,ang_vel,IC9_)
      IMPLICIT NONE
      REAL*8 theta, ang_vel
      INTEGER IR_, IC9_
      CALL GWselvp(IR_, IC9_)
      CALL GWsetogn(IR_, 0, IC9_)
      CALL GWsetpxl(IR_, REAL(theta), REAL(ang_vel), 13)
      END
C
      SUBROUTINE change_amplitude(A,dA,IC1_,IC2_,flag)
      IMPLICIT NONE
      REAL*8 A, dA
      INTEGER k, IC1_, IC2_
      CHARACTER*80 flag
      CALL TBgetkey(k)
      IF(k .EQ. ICHAR('s')) flag = 'stop'
      IF(k .EQ. ICHAR('i')) A = A + dA
      IF(k .EQ. ICHAR('d')) A = A - dA
      IF(k .EQ. ICHAR('i') .OR. k .EQ. ICHAR('d')) THEN
         CALL set_up_windows(A,IC1_,IC2_)
      END IF
      IF(flag .EQ. '') CALL draw_axes(A,IC1_,IC2_)
      END
C
      SUBROUTINE Poincare(theta,ang_vel,IC9_)
*     plot points and allow for change in amplitude of external force
      IMPLICIT NONE
      REAL*8 theta, ang_vel
      INTEGER IR_, IC9_
      CALL GWselvp(IR_, IC9_)
      CALL GWsetogn(IR_, 0, IC9_)
      CALL GWsetpxl(IR_, REAL(theta), REAL(ang_vel), 16)
      END