! motion of a real pendulum with damping and a driving force
! with the Euler-Cromer method
! Program to accompany "Computational Physics" by N. Giordano
! Copyright Prentice Hall 1997
program nonlinear
   option nolet
   library "sgfunc*", "sglib*" 
   dim theta(50000),omega(50000),t(50000)
   call initialize(theta,omega,t,length,dt,damping,force,drive_frequency)
   call calculate(theta,omega,t,length,dt,damping,force,drive_frequency)
   call display(theta,omega,t,length,dt,damping,force,drive_frequency)
end
! initialize variables
! theta() = pendulum angle   omega() = pendulum angular velocity
! damping = damping factor, force = amplitude of drive force
! drive_frequency = angular frequency of drive force
! all angles in radians
sub initialize(theta(),omega(),t(),length,dt,damping,force,drive_frequency)
   input prompt "initial pendulum angle (in radians) -> ": theta(1)
   input prompt "initial angular velocity of pendulum (in radians/s) -> ": omega(1)
   t(1) = 0
   input prompt "length of pendulum (in m) -> ": length
   input prompt "time step -> ": dt
   input prompt "damping constant -> ": damping
   input prompt "amplitude of driving force -> ": force
   input prompt "driving frequency -> ": drive_frequency
end sub
! use the Euler-Cromer method for a damped, nonlinear, driven pendulum
sub calculate(theta(),omega(),t(),length,dt,q,drive_force,drive_frequency)
   i = 0
   g = 9.8
   period = 2 * pi / sqr(g/length)
   do
      i = i + 1
      t(i+1) = t(i) + dt
      omega(i+1) = omega(i) - (g/length) * sin(theta(i)) * dt - q * omega(i) * dt + drive_force * sin(drive_frequency * t(i)) * dt
      theta(i+1) = theta(i) + omega(i+1) * dt
      if theta(i+1) > pi then theta(i+1) = theta(i+1) - 2 * pi ! keep theta in
      if theta(i+1) < -pi then theta(i+1) = theta(i+1) + 2 * pi! range -pi to pi
   loop until t(i+1) >= 60 ! 10 * period
   mat redim omega(i+1),theta(i+1),t(i+1)
end sub
sub display(theta(),omega(),t(),length,dt,damping,force,drive_frequency)
   set background color "white"
   set color "black"
   call setcanvas("white")
   clear
   open #1: screen 0,.9,0,.48
   call settitle("Angular velocity vs. time")
   call sethlabel("Time(s)")
   call setvlabel("Omega")
   call datagraph(t,omega,1,1,"black")
   close #1
   open #1: screen 0,.9,.52,1
   call settitle("Theta vs. time")
   call sethlabel("Time(s)")
   call setvlabel("Theta")
   call datagraph(t,theta,1,1,"black")
   set cursor 3,10
   print "length = "; length
   set cursor 4,10
   print "time step = "; dt
   close #1
end sub