! motion of a simple pendulum   
! theta() = pendulum angle   omega() = pendulum angular velocity t() = time
! length = length of string    dt = time step
! Program to accompany "Computational Physics" by N. Giordano
! Copyright Prentice Hall 1997
program pendulum
   option nolet
   library "sgfunc*","sglib*"
   dim theta(1000), omega(1000), t(1000)
   call initialize(theta,omega,t,length,dt)
   call calculate(theta,omega,t,length,dt)
   call display(theta,omega,t,length,dt)
end
! initialize variables
! theta() = pendulum angle  omega() = pendulum angular velocity
sub initialize(theta(),omega(),t(),length,dt)
   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
end sub
! use the Euler method
sub calculate(theta(),omega(),t(),length,dt)
   i = 0
   g = 9.8
   period = 2 * pi / sqr(g/length)    ! period of pendulum
   do
      i = i + 1
      t(i+1) = t(i) + dt
      omega(i+1) = omega(i) - (g/length) * theta(i) * dt  ! Euler method
      theta(i+1) = theta(i) + omega(i) * dt
   loop until t(i+1) >= 5 * period    ! follow the oscillations for 5 periods
   mat redim omega(i+1),theta(i+1),t(i+1)
end sub
! display theta and omega in separate windows
sub display(theta(),omega(),t(),length,dt)
   set background color "white"
   set color "black"
   call setcanvas("white")
   clear
   open #1: screen 0,.9,0,.48    ! bottom window
   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    ! top window
   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