' gd2lib.bas - jackord@kw.igs.net - revised 28 Feb 08 - Liberty Basic v4.02
' solution of Laplace's equation for concentric square prisms sides L and 2L
'   relaxation calculation with 5120 passes
'   the field is plotted after 20 passes, and whenever the number
'     of passes doubles thereafter
'   the capacitance is calculated using Gauss's law and is expressed as a
'     multiple of that of a cylindrical capacitor with diameters L and 2L
' the relaxation calculation takes a VERY long time to complete

  nomainwin
  button#1.b1, "Plot", [plt], UL, 108, 118, 40, 20
  WindowWidth=274     ' pixel scale 0-256
  WindowHeight=294    ' pixel scale 0-256
  UpperLeftX=100: UpperLeftY=50
  open "Square Prism Capacitor" for graphics_nsb as #1
  #1 "trapclose [quit]"
  #1 "down"
  dim c$(6): dim v(130, 65)
  c$(0)="blue": c$(1)="cyan": c$(2)="red"
  c$(3)="yellow": c$(4)="darkpink": c$(5)="green"
  #1 "color blue ; place 0 0 ; box 257 257"      ' Outer prism
  #1 "color red ; place 64 64 ; box 193 193"     ' Inner prism
  #1 "place 74 86 ; color black"
  #1 "\\Pass    /5120"
  #1 "place 74 150"
  #1 "\\Csqr/Ccir"
  n=64

[waitHere]
  wait

[plt]
  redim v(130, 65)
  for i=n to 2*n: v(i, n)=100: next i                  ' Define V on boundary
  for j=1 to n-1                                       ' Initial linear distn
    for i=j to 2*n: v(i,j)=100*j/n: next i
  next j
  k=1
  for ii=0 to 5121
  #1 "color white ; place 74 86 ; boxfilled 184 106 ; color black"
  #1 "\\Pass "; ii; "/5120"
    for j=1 to n-1                                     ' Project across bdry
      v(j, j+1)=v(j+1, j): v(2*n+1, j)=v(2*n-1, j)
    next j
    if ii=0 or ii=10*k then
      for j=1 to n-1                                   ' Plot potential..
        for i=j to 2*n                                 ' ..in all eight
          plt=int(v(i, j)/7): plt=plt-6*int(plt/6)     ' ..symmetry-related
          #1 "color "; c$(plt)                   ' ..regions
          #1 "line "; i; " "; j; " "; i+1; " "; j
          #1 "line "; 4*n-i; " "; j; " "; 4*n-i+1; " "; j
          #1 "line "; j; " "; i; " "; j+1; " "; i
          #1 "line "; j; " "; 4*n-i; " "; j+1; " "; 4*n-i
          #1 "line "; 4*n-j; " "; i; " "; 4*n-j+1; " "; i
          #1 "line "; 4*n-j; " "; 4*n-i; " "; 4*n-j+1; " "; 4*n-i
          #1 "line "; i; " "; 4*n-j; " "; i+1; " "; 4*n-j
          #1 "line "; 4*n-i; " "; 4*n-j; " "; 4*n-i+1; " "; 4*n-j
        next i
      next j
      j=n/2: s=0                                       ' Calculate C
      for i=j to 2*n                                   ' Gauss's law
        e=v(i, j)-v(i, j-1): s=s+e
      next i
      s=s-e/2: s=s*8/100/16/atn(1): s=s*2*log(2)
      #1 "color white ; place 74 150 ; boxfilled 184 170 ; color black"
      #1 "\\Csq/Ccr "; using("#.###",s)                ' Display C ratio
      k=2*k                                            ' Double plot interval
    end if
    for j=1 to n-1                                     ' The calculation!
      for i=j to 2*n
        v(i, j)=(v(i, j-1)+v(i, j+1)+v(i-1, j)+v(i+1, j))/4
      next i
    next j
    scan
  next ii
goto [waitHere]

[quit]
  close#1
end