PROGRAM poincare
! plot phase diagram and Poincare map for damped, driven pendulum
CALL initial(theta,ang_vel,gamma,A,t,dt,nstep,dA,#1,#2,flag$)
CALL set_up_windows(A,bx,by,#1,#2)
CALL draw_axes(A,#1,#2)
DO
   FOR istep = 1 to nstep
       CALL RK4(theta,ang_vel,gamma,A,t,dt)
       CALL phase_space(theta,ang_vel,#1)
   NEXT istep
   CALL Poincare(theta,ang_vel,bx,by,#2)
   IF key input then CALL change_amplitude(A,dA,#1,#2,flag$)
LOOP UNTIL flag$ = "stop"
END

SUB initial(theta0,ang_vel0,gamma,A,t0,dt,nstep,dA,#1,#2,flag$)
    INPUT prompt "initial angle = ": theta0
    INPUT prompt "initial angular velocity = ": ang_vel0
    LET t0 = 0               ! initial time
    INPUT prompt "damping constant = ": gamma
    INPUT prompt "amplitude of external force = ": A
    ! angular frequency of external force equals two and
    ! hence period of external force equals pi
    LET nstep = 100          ! # iterations between Poincare plot
    LET T = pi               ! period of external force
    LET dt = T/nstep         ! time step
    LET dA = 0.01            ! increase in amplitude of external force
    ! left-side of screen for phase space plot
    OPEN #1: screen 0,0.49,0,1
    ! open second window for Poincare map
    OPEN #2: screen 0.51,1.0,0,1
    LET flag$ = ""
END SUB

SUB set_up_windows(A,bx,by,#1,#2)
    LET vmax = 4*A
    WINDOW #1
    SET WINDOW -pi,pi,-vmax,vmax
    WINDOW #2
    SET WINDOW -pi,pi,-vmax,vmax
    LET scale = 2*min(pi,vmax)
    LET bx = 0.005*scale     ! x-dimension of box
    ASK PIXELS px,py
    IF px > py then
       LET by = (px/py)*bx
    ELSE
       LET by = (py/px)*bx
    END IF
END SUB

SUB draw_axes(A,#1,#2)
    WINDOW #1
    CLEAR
    SET COLOR "black/white"
    LET vmax = 4*A
    PLOT LINES: 0,-vmax;0,vmax
    PLOT LINES: -pi,0;pi,0
    PRINT using "A = #.#####": A
    WINDOW #2
    CLEAR
    SET COLOR "black/white"
    PLOT LINES: -pi,0;pi,0
    PLOT LINES: 0,-vmax; 0,vmax
END SUB

SUB RK4(theta,ang_vel,gamma,A,t,dt)
    ! fourth-order Runge Kutta algorithm
    DECLARE DEF force
    LET k1v = force(theta,ang_vel,gamma,A,t)*dt
    LET k1x = ang_vel*dt
    LET t = t + 0.5*dt
    LET k2v = force(theta+0.5*k1x,ang_vel+0.5*k1v,gamma,A,t)*dt
    LET k2x = (ang_vel + 0.5*k1v)*dt
    LET k3v = force(theta+0.5*k2x,ang_vel+0.5*k2v,gamma,A,t)*dt
    LET k3x = (ang_vel + 0.5*k2v)*dt
    LET t = t + 0.5*dt
    LET k4v = force(theta+k3x,ang_vel+k3v,gamma,A,t)*dt
    LET k4x = (ang_vel + k3v)*dt
    LET ang_vel = ang_vel + (k1v + 2*k2v + 2*k3v + k4v)/6
    LET theta = theta + (k1x + 2*k2x + 2*k3x + k4x)/6
    IF theta > pi then LET theta = theta - 2*pi
    IF theta < -pi then LET theta = theta + 2*pi
END SUB

DEF force(theta,ang_vel,gamma,A,t)
    ! A is amplitude of driving force
    LET force = -gamma*ang_vel - (1 +2*A*cos(2*t))*sin(theta)
END DEF

SUB phase_space(theta,ang_vel,#9)
    WINDOW #9
    SET COLOR "red"
    PLOT theta,ang_vel
END SUB

SUB change_amplitude(A,dA,#1,#2,flag$)
    GET KEY k                ! press any key to clear screen
    IF k = ord("s") then LET flag$ = "stop"
    IF k = ord("i") then LET A = A + dA
    IF k = ord("d") then LET A = A - dA
    IF k = ord("i") or k = ord("d") then
       CALL set_up_windows(A,bx,by,#1,#2)
    END IF
    IF flag$ = "" then CALL draw_axes(A,#1,#2)
END SUB

SUB Poincare(theta,ang_vel,bx,by,#9)
    ! plot points and allow for change in amplitude of external force
    WINDOW #9
    SET COLOR "blue"
    BOX AREA theta-bx,theta+bx,ang_vel-by,ang_vel+by
END SUB