      PROGRAM maxwell
*     implementation of Yee-Visscher algorithm
      IMPLICIT NONE
      REAL*8 bscale, escale, jscale,
     +       B, E, dl, dl_1, dt, fpi, j
      INTEGER IC1_, IC2_, IC3_, IR_, k, mid, n
      LOGICAL Ltmp1_, TBkeyinput
      COMMON /GLBLD/B(0:21,0:21,0:21,3), E(0:21,0:21,0:21,3),
     +             dl, dl_1, dt, fpi, j(0:21,0:21,0:21,3)
      COMMON /GLBLI/mid, n(3)
      DATA IC1_/1/
      DATA IC2_/2/
      DATA IC3_/3/
      CALL GWopen(IR_, 0)
      CALL initial(escale,bscale,jscale)
      CALL screen(IC1_,IC2_,IC3_)
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL current
         CALL newE
         CALL newB
         CALL plotfields(escale,bscale,jscale,IC1_,IC2_,IC3_)
         DO WHILE(.NOT.TBkeyinput())
         END DO
         CALL TBgetkey(k)
         Ltmp1_ = k .NE. ICHAR('s')
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(escale,bscale,jscale)
      IMPLICIT NONE
      REAL*8 pi, escale, bscale, jscale,
     +       B, E, dl, dl_1, dt, fpi, j
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER comp, x, y, z, mid, n
      COMMON /GLBLD/B(0:21,0:21,0:21,3), E(0:21,0:21,0:21,3),
     +             dl, dl_1, dt, fpi, j(0:21,0:21,0:21,3)
      COMMON /GLBLI/mid, n(3)
      dt = 0.03D0
      n(1) = 8
      n(2) = 8
      n(3) = 8
      mid = n(1)/2
      dl = 0.1D0
      dl_1 = 1/dl
      fpi = 4*pi
      escale = dl/(4*pi*dt)
      bscale = escale*dl/dt
      jscale = 1
      DO x = 1, n(1)
         DO y = 1, n(2)
            DO z = 1, n(3)
               DO comp = 1, 3
                  E(x,y,z,comp) = 0
                  B(x,y,z,comp) = 0
               END DO
            END DO
         END DO
      END DO
      END
C
      SUBROUTINE current
*     steady current loop in x-y plane turned on at t = 0 and left on
      IMPLICIT NONE
      REAL*8 B, E, dl, dl_1, dt, fpi, j
      INTEGER mid, n
      COMMON /GLBLD/B(0:21,0:21,0:21,3), E(0:21,0:21,0:21,3),
     +             dl, dl_1, dt, fpi, j(0:21,0:21,0:21,3)
      COMMON /GLBLI/mid, n(3)
      j(mid,mid,mid,2) = 1
      j(mid,mid,mid,1) = -1
      j(mid-1,mid,mid,2) = -1
      j(mid,mid-1,mid,1) = 1
      END
C
      SUBROUTINE newE
*     E defined at the faces
      IMPLICIT NONE
      REAL*8 curlBx, curlBy, curlBz,
     +       B, E, dl, dl_1, dt, fpi, j
      INTEGER x, y, z, mid, n
      COMMON /GLBLD/B(0:21,0:21,0:21,3), E(0:21,0:21,0:21,3),
     +             dl, dl_1, dt, fpi, j(0:21,0:21,0:21,3)
      COMMON /GLBLI/mid, n(3)
      DO x = 1, n(1)
         DO y = 1, n(2)
            DO z = 1, n(3)
               curlBx = B(x,y,z,2)+B(x,y+1,z,3)-B(x,y,z+1,2)-B(x,y,z,3)
               curlBx = curlBx*dl_1
               E(x,y,z,1) = E(x,y,z,1) + dt*(curlBx - fpi*j(x,y,z,1))
               curlBy = B(x,y,z,3)-B(x,y,z,1)+B(x,y,z+1,1)-B(x+1,y,z,3)
               curlBy = curlBy*dl_1
               E(x,y,z,2) = E(x,y,z,2) + dt*(curlBy - fpi*j(x,y,z,2))
               curlBz = B(x,y,z,1)+B(x+1,y,z,2)-B(x,y+1,z,1)-B(x,y,z,2)
               curlBz = curlBz*dl_1
               E(x,y,z,3) = E(x,y,z,3) + dt*(curlBz - fpi*j(x,y,z,3))
            END DO
         END DO
      END DO
      END
C
      SUBROUTINE newB
*     B defined at the edges
      IMPLICIT NONE
      REAL*8 curlEx, curlEy, curlEz,
     +       B, E, dl, dl_1, dt, fpi, j
      INTEGER x, y, z, mid, n
      COMMON /GLBLD/B(0:21,0:21,0:21,3), E(0:21,0:21,0:21,3),
     +             dl, dl_1, dt, fpi, j(0:21,0:21,0:21,3)
      COMMON /GLBLI/mid, n(3)
      DO x = 1, n(1)
         DO y = 1, n(2)
            DO z = 1, n(3)
               curlEx = E(x,y,z,3)-E(x,y,z,2)-E(x,y-1,z,3)+E(x,y,z-1,2)
               curlEx = curlEx*dl_1
               B(x,y,z,1) = B(x,y,z,1) - curlEx*dt
               curlEy = E(x,y,z,1)-E(x,y,z,3)-E(x,y,z-1,1)+E(x-1,y,z,3)
               curlEy = curlEy*dl_1
               B(x,y,z,2) = B(x,y,z,2) - curlEy*dt
               curlEz = E(x,y,z,2)-E(x,y,z,1)-E(x-1,y,z,2)+E(x,y-1,z,1)
               curlEz = curlEz*dl_1
               B(x,y,z,3) = B(x,y,z,3) - curlEz*dt
            END DO
         END DO
      END DO
      END
C
      SUBROUTINE screen(IC1_,IC2_,IC3_)
      IMPLICIT NONE
      REAL*8 B, E, dl, dl_1, dt, fpi, j
      REAL xwin, ywin, GWaspect
      INTEGER IR_, IC1_, IC2_, IC3_, IC4_, mid, n
      COMMON /GLBLD/B(0:21,0:21,0:21,3), E(0:21,0:21,0:21,3),
     +             dl, dl_1, dt, fpi, j(0:21,0:21,0:21,3)
      COMMON /GLBLI/mid, n(3)
      DATA IC4_/4/
      CALL compute_aspect_ratio(REAL(n(1)+1),xwin,ywin)
      CALL GWcolor(IR_, 0, 2)
      CALL GWcolor(IR_, 19, 7)
      CALL GWvport(IR_, 0.0, 0.0, 0.5*GWaspect(1), 0.5)
      CALL GWindow(IR_, 0.0, 0.0, xwin, ywin)
      CALL GWsavevp(IR_, IC1_)
      CALL GWsetogn(IR_, 0, IC1_)
      CALL GWvport(IR_, 0.5*GWaspect(1), 0.5, GWaspect(1), 1.0)
      CALL GWindow(IR_, 0.0, 0.0, xwin, ywin)
      CALL GWsavevp(IR_, IC2_)
      CALL GWsetogn(IR_, 0, IC2_)
      CALL GWvport(IR_, 0.5*GWaspect(1), 0.0, GWaspect(1), 0.5)
      CALL GWindow(IR_, 0.0, 0.0, xwin, ywin)
      CALL GWsavevp(IR_, IC3_)
      CALL GWsetogn(IR_, 0, IC3_)
      CALL GWvport(IR_, 0.0, 0.5, 0.5*GWaspect(1), 1.0)
      CALL GWindow(IR_, 0.0, 24.0, 80.0, 0.0)
      CALL GWsavevp(IR_, IC4_)
      CALL GWsetogn(IR_, 0, IC4_)
      CALL GWputtxt(IR_, 1.0, 2.0, 'Type s to stop')
      CALL GWputtxt(IR_, 1.0, 4.0, 
     +              'Type any other key for next time step')
      END
C
      SUBROUTINE plotfields(escale,bscale,jscale,IC1_,IC2_,IC3_)
      IMPLICIT NONE
      REAL*8 pi, escale, bscale, jscale,
     +       B, E, dl, dl_1, dt, fpi, j
      PARAMETER(pi = 3.141592653589793D0)
      REAL X1, Y1, X2, Y2
      INTEGER IR_, x, y, z, IC1_, IC2_, IC3_, mid, n
      COMMON /GLBLD/B(0:21,0:21,0:21,3), E(0:21,0:21,0:21,3),
     +             dl, dl_1, dt, fpi, j(0:21,0:21,0:21,3)
      COMMON /GLBLI/mid, n(3)
      CALL GWselvp(IR_, IC1_)
      CALL GWsetogn(IR_, 0, IC1_)
      CALL GWerase(IR_, IC1_, -1)
      CALL GWgetwn(IR_, X1, Y1, X2, Y2)
      CALL GWputtxt(IR_, X1+(X2-X1)/10, Y2-(Y2-Y1)/10, 'E(x,y)')
      DO x = 1, n(1)
         DO y = 1, n(2)
            CALL plotarrow(E(x,y,mid,1),x,y,escale,
     +                     0.0D0,0.5D0,0.0D0,pi)
            CALL plotarrow(E(x,y,mid,2),x,y,escale,
     +                     0.5D0,0.0D0,pi/2,3*pi/2)
         END DO
      END DO
      CALL GWselvp(IR_, IC2_)
      CALL GWsetogn(IR_, 0, IC2_)
      CALL GWerase(IR_, IC2_, -1)
      CALL GWgetwn(IR_, X1, Y1, X2, Y2)
      CALL GWputtxt(IR_, X1+(X2-X1)/10, Y2-(Y2-Y1)/10, 'B(x,z)')
      DO x = 1, n(1)
         DO z = 1, n(3)
            CALL plotarrow(B(x,mid,z,1),x,z,bscale,
     +                     0.5D0,0.0D0,0.0D0,pi)
            CALL plotarrow(B(x,mid,z,3),x,z,bscale,
     +                     0.0D0,0.5D0,pi/2,3*pi/2)
         END DO
      END DO
      CALL GWselvp(IR_, IC3_)
      CALL GWsetogn(IR_, 0, IC3_)
      CALL GWerase(IR_, IC3_, -1)
      CALL GWgetwn(IR_, X1, Y1, X2, Y2)
      CALL GWputtxt(IR_, X1+(X2-X1)/10, Y2-(Y2-Y1)/10, 'j(x,y)')
      DO x = 1, n(1)
         DO y = 1, n(2)
            CALL plotarrow(j(x,y,mid,1),x,y,jscale,
     +                     0.0D0,0.5D0,0.0D0,pi)
            CALL plotarrow(j(x,y,mid,2),x,y,jscale,
     +                     0.5D0,0.0D0,pi/2,3*pi/2)
         END DO
      END DO
      END
C
      SUBROUTINE plotarrow(V, x, y, scale, shiftx, shifty, 
     +                     angle1, angle2)
      REAL*8 V, scale, shiftx, shifty, angle1, angle2
      REAL x0(5), y0(5), xx(5), yy(5)
      INTEGER i, x, y
      DATA x0/-0.25, 0.25, 0.12, 0.25, 0.12/, 
     +     y0/ 0.0,  0.0,  0.12, 0.0, -0.12/
      DO i = 1, 5
         xx(i) = x0(i)
         yy(i) = y0(i)
      ENDDO
      IF(V .GT. 0.0D0) THEN
         DO i = 1, 5
            xx(i) = xx(i) * V / scale
            yy(i) = yy(i) * V / scale
            CALL rotate(xx(i), yy(i), REAL(angle1))
            xx(i) = xx(i) + x + shiftx
            yy(i) = yy(i) + y + shifty
         ENDDO
      ELSE IF(V .LT. 0.0D0) THEN
         DOi = 1, 5
            xx(i) = -xx(i) * V / scale
            yy(i) = -yy(i) * V / scale
            CALL rotate(xx(i), yy(i), REAL(angle2))
            xx(i) = xx(i) + x + shiftx
            yy(i) = yy(i) + y + shifty
         ENDDO
      ENDIF
      CALL GWsetpen(IR_, 15, -1, -1, -1)
      CALL GWmove2(IR_, xx(1), yy(1))
      CALL GWline2(IR_, xx(2), yy(2))
      CALL GWline2(IR_, xx(3), yy(3))
      CALL GWmove2(IR_, xx(4), yy(4))
      CALL GWline2(IR_, xx(5), yy(5))
      CALL GWsetpen(IR_, 19, -1, -1, -1)
      END
C
      SUBROUTINE rotate(x, y, angle)
      IMPLICIT NONE
      REAL x, y, angle, tx
      tx = x*cos(angle) - y*sin(angle)
      y = x*sin(angle) + y*cos(angle)
      x = tx
      END
