PROGRAM maxwell
! implementation of Yee-Visscher algorithm
LIBRARY "csgraphics"
PUBLIC E(0 to 21,0 to 21,0 to 21,3),B(0 to 21,0 to 21,0 to 21,3)
PUBLIC j(0 to 21,0 to 21,0 to 21,3)
PUBLIC n(3),dt,dl
PUBLIC mid,fpi,dl_1
CALL initial(escale,bscale,jscale)
CALL screen(#1,#2,#3)
DO
  CALL current(t)
  CALL newE
  CALL newB
  CALL plotfields(escale,bscale,jscale,#1,#2,#3)
  DO
  LOOP until key input
  GET KEY k
LOOP until k = ord("s")
END

SUB initial(escale,bscale,jscale)
   DECLARE PUBLIC E(,,,),B(,,,)
   DECLARE PUBLIC n(),dt
   DECLARE PUBLIC mid,fpi,dl_1
   LET dt = 0.03
   LET n(1) = 8
   LET n(2) = 8
   LET n(3) = 8
   LET mid = n(1)/2
   LET dl = 0.1
   LET dl_1 = 1/dl
   LET fpi = 4*pi
   LET escale = dl/(4*pi*dt)
   LET bscale = escale*dl/dt
   LET jscale = 1
   FOR x = 1 to n(1)
     FOR y = 1 to n(2)
       FOR z = 1 to n(3)
         FOR comp = 1 to 3
           LET E(x,y,z,comp) = 0
           LET B(x,y,z,comp) = 0
         NEXT comp
       NEXT z
     NEXT y
   NEXT x
END SUB

SUB current(t)
   DECLARE PUBLIC j(,,,),n(),mid
   ! steady current loop in x-y plane turned on at t = 0 and left on
   LET j(mid,mid,mid,2) = 1
   LET j(mid,mid,mid,1) = -1
   LET j(mid-1,mid,mid,2) = -1
   LET j(mid,mid-1,mid,1) = 1
END SUB

SUB newE
   ! E defined at the faces
   DECLARE PUBLIC E(,,,),B(,,,),j(,,,)
   DECLARE PUBLIC n(),dt,dl_1,fpi
   FOR x = 1 to n(1)
     FOR y = 1 to n(2)
       FOR z = 1 to n(3)
         LET curlBx = B(x,y,z,2)+B(x,y+1,z,3)-B(x,y,z+1,2)-B(x,y,z,3)
         LET curlBx = curlBx*dl_1
         LET E(x,y,z,1) = E(x,y,z,1) + dt*(curlBx - fpi*j(x,y,z,1))
         LET curlBy = B(x,y,z,3)-B(x,y,z,1)+B(x,y,z+1,1)-B(x+1,y,z,3)
         LET curlBy = curlBy*dl_1
         LET E(x,y,z,2) = E(x,y,z,2) + dt*(curlBy - fpi*j(x,y,z,2))
         LET curlBz = B(x,y,z,1)+B(x+1,y,z,2)-B(x,y+1,z,1)-B(x,y,z,2)
         LET curlBz = curlBz*dl_1
         LET E(x,y,z,3) = E(x,y,z,3) + dt*(curlBz - fpi*j(x,y,z,3))
       NEXT z
     NEXT y
   NEXT x
END SUB

SUB newB
   ! B defined at the edges
   DECLARE PUBLIC E(,,,),B(,,,)
   DECLARE PUBLIC n(),dt,dl_1
   FOR x = 1 to n(1)
     FOR y = 1 to n(2)
       FOR z = 1 to n(3)
         LET curlEx = E(x,y,z,3)-E(x,y,z,2)-E(x,y-1,z,3)+E(x,y,z-1,2)
         LET curlEx = curlEx*dl_1
         LET B(x,y,z,1) = B(x,y,z,1) - curlEx*dt
         LET curlEy = E(x,y,z,1)-E(x,y,z,3)-E(x,y,z-1,1)+E(x-1,y,z,3)
         LET curlEy = curlEy*dl_1
         LET B(x,y,z,2) = B(x,y,z,2) - curlEy*dt
         LET curlEz = E(x,y,z,2)-E(x,y,z,1)-E(x-1,y,z,2)+E(x,y-1,z,1)
         LET curlEz = curlEz*dl_1
         LET B(x,y,z,3) = B(x,y,z,3) - curlEz*dt
       NEXT z
     NEXT y
   NEXT x
END SUB

SUB screen(#1,#2,#3)
   DECLARE PUBLIC n()
   LET L = n(1)
   CALL compute_aspect_ratio(L+1,xwin,ywin)
   SET BACKGROUND COLOR "black"
   OPEN #1: screen 0,.5,0,.5
   SET WINDOW 0,xwin,0,ywin
   SET COLOR "white"
   OPEN #2: screen 0.5,1,0.5,1
   SET WINDOW 0,xwin,0,ywin
   SET COLOR "white"
   OPEN #3: screen 0.5,1,0,0.5
   SET WINDOW 0,xwin,0,ywin
   SET COLOR "white"
   OPEN #4: screen 0,0.5,0.5,1
   SET COLOR "white"
   SET CURSOR 1,1
   PRINT "Type s to stop"
   PRINT
   PRINT "Type any other key for next time step"
END SUB

SUB plotfields(escale,bscale,jscale,#1,#2,#3)
   DECLARE PUBLIC E(,,,),B(,,,),j(,,,)
   DECLARE PUBLIC n(),dt,mid
   WINDOW #1
   CLEAR
   PRINT "E(x,y)"
   FOR x = 1 to n(1)
     FOR y = 1 to n(2)
       CALL plotarrow(E(x,y,mid,1),x,y,escale,0,0.5,0,pi)
       CALL plotarrow(E(x,y,mid,2),x,y,escale,0.5,0,pi/2,3*pi/2)
     NEXT y
   NEXT x
   WINDOW #2
   CLEAR
   PRINT "B(x,z)"
   FOR x = 1 to n(1)
     FOR z = 1 to n(3)
       CALL plotarrow(B(x,mid,z,1),x,z,bscale,0.5,0,0,pi)
       CALL plotarrow(B(x,mid,z,3),x,z,bscale,0,0.5,pi/2,3*pi/2)
     NEXT z
   NEXT x
   WINDOW #3
   CLEAR
   PRINT "j(x,y)"
   FOR x = 1 to n(1)
     FOR y = 1 to n(2)
       CALL plotarrow(j(x,y,mid,1),x,y,jscale,0,0.5,0,pi)
       CALL plotarrow(j(x,y,mid,2),x,y,jscale,0.5,0,pi/2,3*pi/2)
     NEXT y
   NEXT x
END SUB

SUB plotarrow(V,x,y,scale,shiftx,shifty,angle1,angle2)
   IF V > 0 then
      DRAW arrow(V/scale) with rotate(angle1)*shift(x+shiftx,y+shifty)
   ELSE IF V < 0 then
      DRAW arrow(-V/scale) with rotate(angle2)*shift(x+shiftx,y+shifty)
   END IF
END SUB

PICTURE arrow(x)
   SET COLOR "yellow"
   PLOT LINES: -0.25*x,0;0.25*x,0;0.12*x,0.12*x
   PLOT LINES: 0.25*x,0;0.12*x,-0.12*x
   SET COLOR "white"
END PICTURE