*DECK PHYS27
      SUBROUTINE INITAL
C...
C...  ONE-DIMENSIONAL, TIME-DEPENDENT MAXWELL'S EQUATIONS
C...
C...  THE TIME-DEPENDENT MAXWELL'S EQUATIONS ARE (1)
C...
C...        2
C...     DEL E = MU*EPS*E   + MU*SIGMA*E                        (1)
C...                     TT             T
C...
C...        2
C...     DEL H = MU*EPS*H   + MU*SIGMA*H                        (2)
C...                     TT             T
C...
C...  (1)  WYLIE, C. R. AND L. C. BARRETT, ADVANCED ENGINEERING
C...       MATHEMATICS, MCGRAW-HILL BOOK COMPANY, NY, FIFTH EDITION,
C...       1982, PP. 817-820
C...
C...  WHERE
C...
C...     E         ELECTRIC INTENSITY
C...
C...     H         MAGNETIC INTENSITY
C...
C...     T         TIME
C...
C...     DEL       THE LAPLACIAN OPERATOR (A SECOND-ORDER SPATIAL
C...               DIFFERENTIATION OPERATOR)
C...
C...     MU        PERMEABILITY
C...
C...     EPS       PERMITTIVITY
C...
C...     SIGMA     CONDUCTIVITY
C...
C...  AN ALPHABETIC SUBSCRIPT INDICATES A PARTIAL DERIVATIVE, E.G.,
C...  E  AND E   ARE THE FIRST AND SECOND-ORDER PARTIAL DERIVATIVES
C...   T      TT
C...  OF E WITH RESPECT TO T, RESPECTIVELY.
C...
C...  FOR ONE SPATIAL DIMENSION, X, IN CARTESIAN COORDINATES, EQUA-
C...  TIONS (1) AND (2) BECOME
C...
C...     MU*EPS*E   + MU*SIGMA*E  = E                           (3)
C...             TT             T    XX
C...
C...     MU*EPS*H   + MU*SIGMA*H  = H                           (4)
C...             TT             T    XX
C...
C...  CONSIDER EQUATION (3) FOR TWO CASES
C...
C...  (1)  EPS = 0 (GOOD CONDUCTOR)
C...
C...  COMPUTE A SOLUTION TO EQUATION (3) FOR THE INITIAL CONDITION
C...
C...     E(X,0) = 0                                            (5)
C...
C...  AND BOUNDARY CONDITIONS
C...
C...     E(0,T) = 1, E (1,T) = 0                            (6)(7)
C...                  X
C...
C...  USE SUBROUTINE DSS042 FOR THIS SOLUTION.  PLOT E(X,T) VS X FOR
C...  A SERIES OF VALUES OF T.  CHECK THAT YOUR SOLUTION APPROACHES
C...  THE ANALYTICAL SOLUTION FOR LARGE T.  TAKE MU = SIGMA = 1.
C...
C...  (2)  EPS = MU = SIGMA = 1
C...
C...  COMPUTE A SOLUTION TO EQUATION (3) FOR THE INITIAL CONDITIONS
C...
C...     E(X,0) = SIN(PI*X), E (X,0) = 0                    (8)(9)
C...                          T
C...
C...  AND BOUNDARY CONDITIONS
C...
C...     E (0,T) = E (1,T) = 0                            (10)(11)
C...      X         X
C...
C...  USE SUBROUTINES DSS042 AND DSS044 FOR THIS SOLUTION.  PLOT
C...  E(X,T) VS X FOR A SERIES OF VALUES OF T.
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=202)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  U(NX)
     2      /F/ UT(NX)
C...
C...  INITIAL CONDITION(S)
C...
C...     CASE 1
         IF((NORUN.EQ.1).OR.(NORUN.EQ.2))CALL INTAL1
C...
C...     CASE 2
         IF((NORUN.EQ.3).OR.(NORUN.EQ.4))CALL INTAL2
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=202)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  U(NX)
     2      /F/ UT(NX)
C...
C...     CASE 1
         IF((NORUN.EQ.1).OR.(NORUN.EQ.2))CALL DERV1
C...
C...     CASE 2
         IF((NORUN.EQ.3).OR.(NORUN.EQ.4))CALL DERV2
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=202)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  U(NX)
     2      /F/ UT(NX)
C...
C...  PRINT AND PLOT THE SOLUTION
C...
C...     CASE 1
         IF((NORUN.EQ.1).OR.(NORUN.EQ.2))CALL PRINT1(NI,NO)
C...
C...     CASE 2
         IF((NORUN.EQ.3).OR.(NORUN.EQ.4))CALL PRINT2(NI,NO)
      RETURN
      END
      SUBROUTINE INTAL1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=51)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     E(NX)
     2      /F/    ET(NX)
     3      /S/    EX(NX),   EXX(NX)
     4      /C/       EPS,        MU,     SIGMA
     5      /I/         N,        IP
      DOUBLE PRECISION MU
C...
C...  PARAMETERS
      EPS=0.0D0
      MU=1.0D0
      SIGMA=1.0D0
C...
C...  INITIAL CONDITION
      DO 1 I=1,NX
         E(I)=0.0D0
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV1
      IP=0
      RETURN
      END
      SUBROUTINE DERV1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=51)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     E(NX)
     2      /F/    ET(NX)
     3      /S/    EX(NX),   EXX(NX)
     4      /C/       EPS,        MU,     SIGMA
     5      /I/         N,        IP
      DOUBLE PRECISION MU
C...
C...  BOUNDARY CONDITION AT X = 0
       E(1)=1.0D0
      ET(1)=0.0D0
C...
C...  BOUNDARY CONDITION AT X = 1
      EX(NX)=0.0D0
C...
C...  UXX
      IF(NORUN.EQ.1)THEN
         CALL DSS042(0.0D0,1.0D0,NX,E,EX,EXX,1,2)
      ELSE IF(NORUN.EQ.2)THEN
         CALL DSS044(0.0D0,1.0D0,NX,E,EX,EXX,1,2)
      END IF
C...
C...  PDE
      DO 1 I=2,NX
         ET(I)=EXX(I)
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT1(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=51)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     E(NX)
     2      /F/    ET(NX)
     3      /S/    EX(NX),   EXX(NX)
     4      /C/       EPS,        MU,     SIGMA
     5      /I/         N,        IP
      DOUBLE PRECISION MU
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      COMMON/P1/   XP(NX),  EP(6,NX)
C...
C...  PRINT A HEADING FOR THE SOLUTION
C...
C...  MONITOR THE OUTPUT
      WRITE(*,*)NORUN,T
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IP=IP+1
      IF(IP.EQ.1)WRITE(NO,1)
1     FORMAT(9X,'T',/,
     1 3X,'   E(0,T)',3X,' E(0.2,T)',3X,' E(0.4,T)',
     2 3X,' E(0.6,T)',3X,' E(0.8,T)',3X,' E(1.0,T)',/)
C...
C...  PRINT THE NUMERICAL SOLUTION
      WRITE(NO,2)T,(E(I),I=1,NX,10)
2     FORMAT(F9.2,/,6F12.5,/)
C...
C...  STORE THE SOLUTION FOR PLOTTING
      DO 3 I=1,NX
         XP(I)=DFLOAT(I-1)/DFLOAT(NX-1)
         EP(IP,I)=E(I)
3     CONTINUE
      CALL PLOT1
C...
C...  PLOT THE SOLUTION
      IF(IP.LT.6)RETURN
      CALL SPLOTS(IP,NX,XP,EP)
      WRITE(NO,4)
4     FORMAT(1H ,/,' E(X,T) VS X, T = 0, 0.2,..., 1.0')
      RETURN
      END
      SUBROUTINE PLOT1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=51)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     E(NX)
     2      /F/    ET(NX)
     3      /S/    EX(NX),   EXX(NX)
     4      /C/       EPS,        MU,     SIGMA
     5      /I/         N,        IP
      DOUBLE PRECISION MU
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      COMMON/P1/   XP(NX),  EP(6,NX)
C...
C...  WRITE A FILE FOR TOP DRAWER PLOTTING
      IF(IP.EQ.1)THEN
         OPEN(2,FILE='T.TOP',STATUS='NEW')
         WRITE(2,12)
12       FORMAT(' SET LIMITS X FROM 0 TO 1 Y FROM 0 TO 1',/,
     1          ' SET FONT DUPLEX')
      END IF
13    FORMAT(' JOIN 1')
C...
C...  WRITE THE NUMERICAL SOLUTION AND PLOT LABELS
      IF(IP.EQ.6)THEN
         WRITE(2,1)
1        FORMAT(' SET SYMBOL 0')
         WRITE(2,11)(XP(I),EP(1,I),I=1,NX)
11       FORMAT(2F10.4)
         WRITE(2,13)
         WRITE(2,2)
2        FORMAT(' SET SYMBOL 2')
         WRITE(2,11)(XP(I),EP(2,I),I=1,NX)
         WRITE(2,13)
         WRITE(2,3)
3        FORMAT(' SET SYMBOL 4')
         WRITE(2,11)(XP(I),EP(3,I),I=1,NX)
         WRITE(2,13)
         WRITE(2,4)
4        FORMAT(' SET SYMBOL 6')
         WRITE(2,11)(XP(I),EP(4,I),I=1,NX)
         WRITE(2,13)
         WRITE(2,5)
5        FORMAT(' SET SYMBOL 8')
         WRITE(2,11)(XP(I),EP(5,I),I=1,NX)
         WRITE(2,13)
         WRITE(2,6)
6        FORMAT(' SET SYMBOL 1')
         WRITE(2,11)(XP(I),EP(6,I),I=1,NX)
         WRITE(2,13)
      WRITE(2,14)NORUN
14    FORMAT(
     +' SET FONT DUPLEX'
     +,/,' TITLE 4.5 9.5',
     +   ' "Figure ',I2,'a: Maxwell Equation, Eps = 0"'
     +,/,' TITLE LEFT "E(x,t)"'
     +,/,' TITLE BOTTOM "x"'
     +,/,' TITLE 6.25 0.5 "E(x,t) vs  x"'
     +,/,' TITLE 4.25 0.10 "t = 0, 0.2,..., 1, NX = 51',
     +   ' Uniform Grid"')
C...
C...  NEXT PLOT
      WRITE(2,15)
15    FORMAT('NEW FRAME')
      END IF
      RETURN
      END
      SUBROUTINE INTAL2
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=101)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    E1(NX),    E2(NX)
     2      /F/   E1T(NX),   E2T(NX)
     3      /S/    EX(NX),   EXX(NX)
     4      /C/       EPS,        MU,     SIGMA
     5      /I/         N,        IP
      DOUBLE PRECISION MU
C...
C...  SET THE PARAMETERS
      EPS=1.0D0
      MU=1.0D0
      SIGMA=1.0D0
C...
C...  INITIAL CONDITION
      PI=3.1415927D0
      DO 1 I=1,NX
         X=DFLOAT(I-1)/DFLOAT(NX-1)
         E1(I)=DSIN(PI*X)
1     CONTINUE
C...
C...  INITIAL CONDITION
      DO 2 I=1,NX
         E2(I)=0.0D0
2     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV2
      IP=0
      RETURN
      END
      SUBROUTINE DERV2
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=101)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    E1(NX),    E2(NX)
     2      /F/   E1T(NX),   E2T(NX)
     3      /S/    EX(NX),   EXX(NX)
     4      /C/       EPS,        MU,     SIGMA
     5      /I/         N,        IP
      DOUBLE PRECISION MU
C...
C...  BOUNDARY CONDITION AT X = 0
      EX(1)=0.0D0
C...
C...  BOUNDARY CONDITION AT X = 1
      EX(NX)=0.0D0
C...
C...  UXX
      IF(NORUN.EQ.3)THEN
         CALL DSS042(0.0D0,1.0D0,NX,E1,EX,EXX,2,2)
      ELSE IF(NORUN.EQ.4)THEN
         CALL DSS044(0.0D0,1.0D0,NX,E1,EX,EXX,2,2)
      END IF
C...
C...  PDES
      DO 1 I=1,NX
         E1T(I)=E2(I)
         E2T(I)=EXX(I)-E2(I)
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT2(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=101)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    E1(NX),    E2(NX)
     2      /F/   E1T(NX),   E2T(NX)
     3      /S/    EX(NX),   EXX(NX)
     4      /C/       EPS,        MU,     SIGMA
     5      /I/         N,        IP
      DOUBLE PRECISION MU
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      COMMON/P2/   XP(NX),  EP(6,NX)
C...
C...  MONITOR THE OUTPUT
      WRITE(*,*)NORUN,T
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IP=IP+1
      IF(IP.EQ.1)WRITE(NO,1)
1     FORMAT(9X,'T',/,
     1 3X,'   E(0,T)',3X,' E(0.2,T)',3X,' E(0.4,T)',
     2 3X,' E(0.6,T)',3X,' E(0.8,T)',3X,' E(1.0,T)',/)
C...
C...  PRINT THE NUMERICAL SOLUTION
      WRITE(NO,2)T,(E1(I),I=1,NX,20)
2     FORMAT(F9.2,/,6F12.5,/)
C...
C...  STORE THE SOLUTION FOR PLOTTING
      DO 3 I=1,NX
         XP(I)=DFLOAT(I-1)/DFLOAT(NX-1)
         EP(IP,I)=E1(I)
3     CONTINUE
      CALL PLOT2
C...
C...  PLOT THE SOLUTION
      IF(IP.LT.6)RETURN
      CALL SPLOTS(IP,NX,XP,EP)
      WRITE(NO,4)
4     FORMAT(1H ,/,' E(X,T) VS X, T = 0, 0.2,..., 1.0')
      RETURN
      END
      SUBROUTINE PLOT2
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=101)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    E1(NX),    E2(NX)
     2      /F/   E1T(NX),   E2T(NX)
     3      /S/    EX(NX),   EXX(NX)
     4      /C/       EPS,        MU,     SIGMA
     5      /I/         N,        IP
      DOUBLE PRECISION MU
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      COMMON/P2/   XP(NX),  EP(6,NX)
C...
C...  WRITE A FILE FOR TOP DRAWER PLOTTING
      IF(IP.EQ.1)THEN
         OPEN(2,FILE='T.TOP',STATUS='NEW')
         WRITE(2,12)
12       FORMAT(' SET LIMITS X FROM 0 TO 1 Y FROM 0 TO 1',/,
     1          ' SET FONT DUPLEX')
      END IF
13    FORMAT(' JOIN 1')
C...
C...  WRITE THE NUMERICAL SOLUTION AND PLOT LABELS
      IF(IP.EQ.6)THEN
         WRITE(2,1)
1        FORMAT(' SET SYMBOL 0')
         WRITE(2,11)(XP(I),EP(1,I),I=1,NX)
11       FORMAT(2F10.4)
         WRITE(2,13)
         WRITE(2,2)
2        FORMAT(' SET SYMBOL 2')
         WRITE(2,11)(XP(I),EP(2,I),I=1,NX)
         WRITE(2,13)
         WRITE(2,3)
3        FORMAT(' SET SYMBOL 4')
         WRITE(2,11)(XP(I),EP(3,I),I=1,NX)
         WRITE(2,13)
         WRITE(2,4)
4        FORMAT(' SET SYMBOL 6')
         WRITE(2,11)(XP(I),EP(4,I),I=1,NX)
         WRITE(2,13)
         WRITE(2,5)
5        FORMAT(' SET SYMBOL 8')
         WRITE(2,11)(XP(I),EP(5,I),I=1,NX)
         WRITE(2,13)
         WRITE(2,6)
6        FORMAT(' SET SYMBOL 1')
         WRITE(2,11)(XP(I),EP(6,I),I=1,NX)
         WRITE(2,13)
      WRITE(2,14)NORUN
14    FORMAT(
     +' SET FONT DUPLEX'
     +,/,' TITLE 4.5 9.5',
     +   ' "Figure ',I2,'a: Maxwell Equation, Eps = 1"'
     +,/,' TITLE LEFT "E(x,t)"'
     +,/,' TITLE BOTTOM "x"'
     +,/,' TITLE 6.25 0.5 "E(x,t) vs  x"'
     +,/,' TITLE 4.25 0.1 "t = 0, 0.2,..., 1, NX = 101',
     +   ' Uniform Grid"')
C...
C...  NEXT PLOT
      WRITE(2,15)
15    FORMAT('NEW FRAME')
      END IF
      RETURN
      END
