      PROGRAM oscillators
*     simulate coupled linear oscillators in one dimension
      IMPLICIT NONE
      REAL*8 dt, ke, t, u(0:21), v(0:21), usave(0:21)
      INTEGER IR_, L, i, loop
      PARAMETER(loop = 4)
      LOGICAL TBkeyinput
      CALL GWopen(IR_, 0)
      CALL initial(L,u,v,t,dt,usave)
      CALL GWsleep2(IR_, 0)
      DO WHILE(.NOT.TBkeyinput())
         DO i = 1, loop
            CALL update(L,u,v,ke,t,dt,usave)
         END DO
         CALL animate(L,u,ke,t,usave)
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(L,u,v,t,dt,usave)
      IMPLICIT NONE
      REAL*8 u(0:21), v(0:21), t, dt, usave(0:21), DATA_(20)
      INTEGER IR_, j, L, Itmp1_
      INTEGER DP_
      DATA DP_/1/, DATA_/0.0D0,0.5D0,18*0.0D0/
      t = 0
      dt = 0.025D0
      L = 10                                ! number of particles
      CALL GWindow(IR_, -1.0, -4.0, REAL(L+1), 4.0)
      CALL GWline(IR_, 0.5, 0.0, L+0.5, 0.0)
      CALL GWanchor(IR_, 1)
      CALL GWsetogn(IR_, 0, 1)
      DO j = 1, L
         u(j) = DATA_(DP_)
         DP_ = DP_ + 1
         CALL GWsrect(IR_, REAL(j), REAL(u(j)), 
     +                     REAL(j)+0.2, REAL(u(j))+0.2, 13)
      END DO
      DO j = 1, L
         v(j) = DATA_(DP_)                  ! initial velocities
         DP_ = DP_ + 1
      END DO
      u(0) = 0                              ! fixed wall boundary conditions
      u(L+1) = 0
      DO Itmp1_ = 0, 21
         usave(Itmp1_) = u(Itmp1_)
      END DO
      END
C
      SUBROUTINE update(L,u,v,ke,t,dt,usave)
*     Euler-Richardson algorithm
      IMPLICIT NONE
      REAL*8 a(20), amid(20), umid(0:21), vmid(20), u(0:21), v(0:21),
     +       ke, t, dt, usave(0:21)
      INTEGER j, L
      ke = 0
*     K/M equal to unity
      DO j = 1, L
         usave(j) = u(j)
         a(j) = -2*u(j) + u(j+1) + u(j-1)
         umid(j) = u(j) + 0.5D0*v(j)*dt
         vmid(j) = v(j) + 0.5D0*a(j)*dt
      END DO
      umid(0) = 0
      umid(L+1) = 0
      DO j = 1, L
         amid(j) = -2*umid(j) + umid(j+1) + umid(j-1)
         u(j) = u(j) + vmid(j)*dt
         v(j) = v(j) + amid(j)*dt
         ke = ke + v(j)*v(j)
      END DO
      t = t + dt
      END
C
      SUBROUTINE animate(L,u,ke,t,usave)
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 E, pe, u(0:21), ke, t, usave(0:21)
      INTEGER IR_, j, L
      pe = (u(1) - u(0))**2                 ! interaction with left spring
      CALL GWsetogn(IR_, 0, -1)
      DO j = 1, L
*        compute potential energy
         pe = pe + (u(j+1) - u(j))**2
*        transverse oscillation
         CALL GWsrect(IR_, REAL(j), REAL(usave(j)), 
     +                     REAL(j)+0.2, REAL(usave(j))+0.2, 13)
      END DO
      CALL GWflush(IR_, -1)
      CALL GWsetogn(IR_, 0, -2)
      WRITE(Stmp1_, '(A,F6.2)') 't = ', t
      CALL GWputtxt(IR_, 0.0, 3.7, Stmp1_)
      E = 0.5D0*(ke + pe)
      WRITE(Stmp1_, '(A,F6.4)') 'E = ', E
      CALL GWputtxt(IR_, 0.0, 3.4, Stmp1_)
      CALL GWflush(IR_, -2)
      END
