! simulate waves on a string
! Program to accompany "Computational Physics" by N. Giordano
! Copyright Prentice Hall 1997
program waves
   option NOLET
   dim x(0 to 1000,3)
   call init(x,dt,dx,c,nmax)
   call first_display(x,dx,nmax)
   i = 0
   do
      call propagate(x,dt,dx,c,nmax)
      call display(x,dx,nmax)
      if key input then 
         get key z
         if chr$(z) = " " then exit do
      end if
   loop
end
!   y(,) is a two dimensional array containing the displacement at
!      the current time step y(i,2), the previous time step y(i,1)
!      the next time step corresponds to y(i,3)
!   dt = time step     dx = spatial grid size   c = wave speed
!   nmax = number of spatial units 
!   i runs from 0 to nmax, with the fixed ends at i=0 and i=nmax
sub propagate(y(,),dt,dx,c,nmax)
   r2 = (dt * c / dx)^2
   for i = 1 to nmax-1
      y(i,3) = r2 * (y(i+1,2) + y(i-1,2)) + 2*(1.0 - r2)*y(i,2) - y(i,1)
   next i
   for i = 1 to nmax-1
      y(i,1) = y(i,2)
      y(i,2) = y(i,3)
   next i
end sub
! display string profile
sub display(y(,),dx,nmax)
   clear
   set color "red"
   plot 0,0;dx*nmax,0     ! plot axes
   set color "black"
   plot 0,y(0,2);
   for i = 1 to nmax-1    ! plot each interior grid point
      plot i * dx,y(i,2);
   next i 	
   plot nmax*dx,y(nmax,2)
end sub
! set up for first plot
sub first_display(y(,),dx,nmax)
   set background color "white"
   clear
   set window 0,nmax*dx,-1.5,1.5
   call display(y,dx,nmax)
end sub
! initialize variables   y(i,n) = string position at site i
! n = 3 = new position
! n = 2 = previous position
! n = 1 = position two time steps ago
! dt = time step   dx = spatial step   c = wave speed  
! nmax = number of grid steps
sub init(y(,),dt,dx,c,nmax)
   input prompt "number of spatial grid points (try 100) -> ": nmax 
   for i = 1 to nmax   ! initial all to zero
      y(i,1) = 0
      y(i,2) = 0
      y(i,3) = 0
   next i
   dx = .01
   c = 300
   dt = dx / c         ! choose so that r = 1 for stability
   for i = 2 to nmax - 1
      y(i,2) = exp(-1000 * (i * dx - 0.3 * nmax * dx)^2)  ! initial profile
      y(i,1) = exp(-1000 * (i * dx - 0.3 * nmax * dx)^2)
   next i
end sub