REM gd2bbc.bas - jackord@kw.igs.net - revised 27 Feb 08 - BBC BASIC (RTR)
      REM solution of Laplace's equation for concentric square prisms sides L and 2L
      REM   relaxation calculation with 5120 passes
      REM   the field is plotted after 20 passes, and whenever the number
      REM     of passes doubles thereafter
      REM   the capacitance is calculated using Gauss's law and is expressed as a
      REM     multiple of that of a cylindrical capacitor with diameters L and 2L
      REM the relaxation calculation takes a VERY long time to complete
      
      *FLOAT 64
      INSTALL @lib$+"WINLIB5"
      
      button% = FN_button("Plot", 108, 118, 40, 20, FN_setproc(PROCplt), 0)
      WindowWidth=256     : REM pixel scale 0-256
      WindowHeight=256    : REM pixel scale 0-256
      VDU 23,22,WindowWidth;WindowHeight;8,15,16,128
      SYS "SetWindowText", @hwnd%, "Square Prism Capacitor"
      
      ON CLOSE PROCundim(v()) : QUIT
      DIM c%(6): PROCdim2d(v(), 8, 130, 65)
      c%(0)=12: c%(1)=14: c%(2)=9
      c%(3)=11: c%(4)=5: c%(5)=10
      GCOL 12 : RECTANGLE 0, 0, 510, 510 : REM Outer prism
      GCOL 9 : RECTANGLE 130, 128, 254, 256 : REM Inner prism
      GCOL 0 : VDU 5 : MOVE 148, 212 : PRINT "Pass    /5120";
      MOVE 148, 340 : PRINT "Csqr/CCir";
      n=64
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      
      DEF PROCplt
      FOR i=n TO 2*n: v(i, n)=100: NEXT i : REM Define V on boundary
      FOR j=1 TO n-1                      : REM 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
        GCOL 15 : RECTANGLE FILL 148, 172, 220, 40
        GCOL 0 : MOVE 148, 212 : PRINT "Pass ";ii; "/5120";
        FOR j=1 TO n-1                    : REM 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                  : REM Plot potential..
            FOR i=j TO 2*n                : REM ..in all eight
              plt=INT(v(i, j)/7): plt=plt-6*INT(plt/6) : REM ..symmetry-related
              GCOL c%(plt)                : REM ..regions
              LINE i*2, j*2, i*2+2, j*2
              LINE 8*n-i*2, j*2, 8*n-2*i+2, j*2
              LINE j*2, i*2, j*2+2, i*2
              LINE j*2, 8*n-i*2, j*2+2, 8*n-i*2
              LINE 8*n-j*2, i*2, 8*n-j*2+2, i*2
              LINE 8*n-j*2, 8*n-i*2, 8*n-j*2+2, 8*n-i*2
              LINE i*2, 8*n-j*2, i*2+2, 8*n-j*2
              LINE 8*n-i*2, 8*n-j*2, 8*n-i*2+2, 8*n-j*2
            NEXT i
          NEXT j
          j=n/2: s=0                      : REM Calculate C
          FOR i=j TO 2*n                  : REM 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)
          GCOL 15 : RECTANGLE FILL 148, 300, 220, 40
          GCOL 0 : MOVE 148, 340
          PRINT "Csq/Ccr "; FNusing("#.###",s); : REM Display C ratio
          k=2*k                           : REM Double PLOT INTerval
        ENDIF
        FOR j=1 TO n-1                    : REM 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
      NEXT ii
      ENDPROC
      
      DEF PROCdim2d(RETURN A(),S%,D1%,D2%)
      LOCAL A%
      SYS "GlobalAlloc", 64, 9+S%*(D1%+1)*(D2%+1) TO A%
      ?A%=2 : A%!1=D1%+1 : A%!5=D2%+1
      !^A() = A%
      ENDPROC
      
      DEF PROCundim(RETURN A())
      SYS "GlobalFree", !^A()
      !^A() = 0
      ENDPROC
      
      REM FNUSING.BBC Version 1.3
      REM Richard Russell 26-Jun-2005
      ;
      DEF FNusing(F$,N):LOCAL @%,I%,F%,L%,P%,R%,S$
      REPEAT I% += 1
        L%=I%
        IF MID$(F$,I%,3)="**$" I% += 3 : F% OR=3
        IF MID$(F$,I%,3)="**£" I% += 3 : F% OR=5
        IF MID$(F$,I%,2)="**"  I% += 2 : F% OR=1
        IF MID$(F$,I%,2)="$$"  I% += 2 : F% OR=2
        IF MID$(F$,I%,2)="££"  I% += 2 : F% OR=4
        IF MID$(F$,I%,1)="+"   I% += 1 : F% OR=8
        WHILE MID$(F$,I%,1)="#" I% += 1
          IF MID$(F$,I%,1)=","  I% += 1 : F% OR=16
        ENDWHILE
        IF MID$(F$,I%,1)="."   I% += 1 : F% OR=32
        P%=I%
        WHILE MID$(F$,I%,1)="#" I% += 1 : ENDWHILE
        IF MID$(F$,I%,4)="^^^^" I% += 4 : F% OR=128
        R%=I%
        IF MID$(F$,I%,1)="+" I% += 1 : F% OR=64
        IF MID$(F$,I%,1)="-" MID$(F$,I%,1)=" " : I% += 1 : F% OR=64
      UNTIL (I%>=LENF$) OR (L%<>R%)
      :
      IF F% AND 128 @%=&1010000+(R%-P%-3)*256 ELSE @%=&1020000+(R%-P%)*256
      S$=STR$(ABSN)
      IF F% AND 16 I%=LENS$-R%+P%+(R%<>P%)-2 : WHILE I%>1 S$=LEFT$(S$,I%-1)+","+MID$(S$,I%) : I%-=3 : ENDWHILE
      IF F% AND 32 IF (P%-L%)=1 IF ASC(S$)=48 S$=MID$(S$,2)
      IF F% AND 2 S$="$"+S$
      IF F% AND 4 S$="£"+S$
      IF F% AND 8 IF N>=0 S$="+"+S$
      IF (F% AND 64)=0 IF N<0 S$="-"+S$
      IF F% AND 1 MID$(F$,L%,R%-L%)=STRING$(R%-L%,"*") ELSE MID$(F$,L%,R%-L%)=STRING$(R%-L%," ")
      I%=LEN(S$)
      IF I%<=(R%-L%) MID$(F$,R%-I%,I%)=S$ ELSE MID$(F$,L%,R%-L%)=STRING$(R%-L%,"?")
      IF F% AND 64 IF N<0 MID$(F$,R%,1)="-"
      =F$