! calculate the electric potential and field lines
! between two parallel plates
! use the Jacobi method
! Program to accompany "Computational Physics" by N. Giordano
! Copyright Prentice Hall 1997
program laplace
   option nolet
   library "sgfunc*", "sglib*" 
   dim v(0,0)    ! store the potential in this array
   call initialize(v,max)
   call calculate(v,max,n)     ! iterate with relaxation method
   call display(v(,),max)    ! plot final result
end
! initialize variables
! v() = potential   i,j = position   max determines number of grid steps
sub initialize(v(,),max)
   print "number of grid steps per unit distance"
   input prompt "(x and y run from -1 to 1) -> ": max
   mat redim v(-max to max,-max to max)  ! resize array
   for i = -max to max                   ! zero out array
      for j = -max to max
         v(i,j) = 0
      next j
   next i
   for j = -max to max                   ! set up boundary conditions
      v(-max,j) = -1                     ! along edges
      v(max,j) = 1
   next j
   for i = -max to max
      v(i,-max) = i / max 
      v(i,max) = i / max
   next i
end sub
! v() is the potential, max -> determines number of grid steps
! n = number of iterations of relaxation algorithm
sub calculate(v(,),max,n)
   dim v_tmp(0,0)
   mat redim v_tmp(-max to max,-max to max)
   for i = -max to max    ! initialize temporary array 
      for j = -max to max       ! not absolutely necessary, but a good
         v_tmp(i,j) = v(i,j)    ! good habit to get into
      next j
   next i
   i = 0
   do                     ! main loop - the real work is done here
      i = i + 1
      call update(v,v_tmp,max,diff,n)
      call update(v_tmp,v,max,diff,n)
   loop until i > 10 and diff < 1e-5
end sub
! use Jacobi relaxation algorithm to calculate new values of the potential
! old values are in v1, put new values into v2
! diff is average change in potential during the current iteration
sub update(v1(,),v2(,),max,diff,n)
   diff = 0
   for i = -max+1 to max-1
      for j = -max+1 to max-1
         tmp = v1(i,j)
         v2(i,j) = (v1(i-1,j) + v1(i+1,j) + v1(i,j+1) + v1(i,j-1)) / 4
         diff = diff + abs(tmp - v2(i,j))
      next j
   next i
   diff = diff / (2*max+1)^2
end sub
! display final results
sub display(v(,),max)
   dim x(0),y(0),ex(0),ey(0)
   set background color "white"
   set color "black"
   clear
   call settitle("Equipotential lines")
   call sethlabel("X")
   call setvlabel("Y")
   call topograph(v,-1,1,-1,1,"black")  ! equipotential lines
   get key z
   max2 = (2*max+1)^2                   ! prepare for vector plot
   mat redim x(max2),y(max2),ex(max2),ey(max2)
   k = 0
   for i = -max+1 to max-1
      for j = -max+1 to max-1
         k = k + 1
         x(k) = i/max
         y(k) = j/max
         ex(k) = -(v(i+1,j) - v(i-1,j)) * max   ! compute electric field
         ey(k) = -(v(i,j+1) - v(i,j-1)) * max   ! components
      next j
   next i
   for i = -max to max step 2*max
      k = k + 1
      x(k) = i/max
      y(k) = i/max
      ex(k) = 0
      ey(k) = 0
   next i
   mat redim x(k),y(k),ex(k),ey(k)
   call settitle("Electric field between two metal plates")
   call sethlabel("x")
   call setvlabel("y")
   call vectorgraph(x,y,ex,ey,.1,3,"black")  ! make vector plot
end sub