      SUBROUTINE INITAL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NR=51)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     C(NR)
     2      /F/    CT(NR)
     3      /S/    CR(NR),   CRR(NR)
     4      /R/        R1,     R(NR),       DR,         PI,      CMAX,
     5                 MU,     SIGMA,       WI
     6      /I/        IP
      DOUBLE PRECISION MU
C...
C...  RADIAL GRID
C...
C...     CONDUCTOR RADIUS
         R1=1.0D0
C...
C...     GRID INCREMENT
         DR=R1/DFLOAT(NR-1)
C...
C...     GRID POINTS
         DO 1 I=1,NR
            R(I)=DFLOAT(I-1)*DR
1        CONTINUE
C...
C...  PI
         PI=4.0D0*DATAN(1.0D0)
C...
C...  MAXIMUM CONDUCTOR CURRENT
         CMAX=1.0D+03
C...
C...  PERMEABILITY
         MU=1.0D0
C...
C...  CONDUCTIVITY
         SIGMA=1.0D+01
C...
C...  CURRENT FREQUENCY
         WI=1.0D0
C...
C...  INITIAL CONDITION
         DO 2 I=1,NR
            C(I)=0.0D0
2        CONTINUE
C...
C...  INITIALIZE COUNTER FOR PRINTED AND PLOTTED SOLUTION
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NR=51)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     C(NR)
     2      /F/    CT(NR)
     3      /S/    CR(NR),   CRR(NR)
     4      /R/        R1,     R(NR),       DR,         PI,      CMAX,
     5                 MU,     SIGMA,       WI
     6      /I/        IP
      DOUBLE PRECISION MU
C...
C...  DERIVATIVE J
C...              R
      CALL DSS004(0.0D0,R1,NR,C,CR)
C...
C...  BOUNDARY CONDITION AT R = 0
      NL=2
      CR(1)=0.0D0
C...
C...  BOUNDARY CONDITION AT R = R1
      NU=2
      CR(NR)=(MU*SIGMA)/(2.0D0*PI*R1)*CMAX*WI*DCOS(WI*T)
C...
C...  DERIVATIVE J
C...              RR
      CALL DSS044(0.0D0,R1,NR,C,CR,CRR,NL,NU)
C...
C...  PDE
      DO 1 I=1,NR
         IF(I.EQ.1)THEN
            CT(I)=(1.0D0/(MU*SIGMA))*2.0D0*CRR(I)
         ELSE
            CT(I)=(1.0D0/(MU*SIGMA))*(CRR(I)+1.0D0/R(I)*CR(I))
         END IF
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NR=51)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     C(NR)
     2      /F/    CT(NR)
     3      /S/    CR(NR),   CRR(NR)
     4      /R/        R1,     R(NR),       DR,         PI,      CMAX,
     5                 MU,     SIGMA,       WI
     6      /I/        IP
      DOUBLE PRECISION MU
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IP=IP+1
      IF(IP.EQ.1)THEN
         WRITE(NO,1)
1        FORMAT(7X,'T'
     1   ,1X,'J(0.0*r1,t)',1X,'J(0.2*r1,t)',1X,'J(0.4*r1,t)',
     2  /,9X,'J(0.6*r1,t)',1X,'J(0.8*r1,t)',1X,'J(1.0*r1,t)',
     3  /,9X,'       itot',1X,'       i(t)',1X,'       diff',
     4  /)
         END IF
C...
C...  INTEGRAL OF THE CURRENT DENSITY BY SIMPSON'S RULE
      CTOT=0.
      DO 2 I=2,NR-1,2
         CTOT=CTOT+(R(I-1)*C(I-1)+4.0D0*R(I)*C(I)+R(I+1)*C(I+1))
2     CONTINUE
      CTOT=2.0D0*PI*DR/3.0D0*CTOT
C...
C...  CONDUCTOR CURRENT AND DIFFERENCE OF CURRENTS
      CC=CMAX*DSIN(WI*T)
      DIFF=CTOT-CC
C...
C...  PRINT THE NUMERICAL SOLUTION
      NRI=10
      WRITE(NO,3)T,(C(I),I=1,NR,NRI),CTOT,CC,DIFF
3     FORMAT(F8.3,3F12.1,/,8X,3F12.1,/,8X,3F12.1,/)
C...
C...  STORE AND PLOT THE SOLUTION
      CALL PLOT
      RETURN
      END
      SUBROUTINE PLOT
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER(NR=51)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/     C(NR)
     2      /F/    CT(NR)
     3      /S/    CR(NR),   CRR(NR)
     4      /R/        R1,     R(NR),       DR,         PI,      CMAX,
     5                 MU,     SIGMA,       WI
     6      /I/        IP
      DOUBLE PRECISION MU
C...
C...  INITIALIZE PLOT NUMBER
      DATA NCASE/0/
C...
C...  OPEN FILE FOR TOP DRAWER PLOTTING
      IF(IP.EQ.1)THEN
         OPEN(4,FILE='T.TOP',STATUS='NEW')
C...
C...     WRITE TOP DRAW FILE FOR PLOTTING CURRENT DENSITY
         WRITE(4,17)
17       FORMAT(' SET LIMITS X FROM 0 TO 1 Y FROM -600 TO 600',/,
     1          ' SET FONT DUPLEX')
C...
C...     SET THE CASE NUMBER FOR THE PLOT LABEL
         NCASE=NCASE+1
      END IF
C...
C...  WRITE THE SOLUTION IN FILE T.TOP FOR PLOTTING
      WRITE(4,14)(R(I),C(I),I=1,NR)
14    FORMAT(2F10.3)
      WRITE(4,16)
16    FORMAT(' JOIN 1')
C...
C...  AT THE END OF THE SOLUTION, LABEL THE PLOT
      IF(IP.EQ.5)THEN
         WRITE(4,18)NCASE
18       FORMAT(
     +   ' Title 4.5 9.5 "Figure ',I2,
     +   ': Current Density - J(r,t)"'
     +   ,/,' TITLE LEFT "J(r,t) (amps/cm**2)"'
     +   ,/,' TITLE BOTTOM "r (cm)"'
     +   ,/,'TITLE 3.5 0.40 '
     +      '"(mu)(sigma) = 10.0 sec/cm**2, wi = 1.0 1/sec"')
C...
C...     NEW PLOT
         WRITE(4,19)
19       FORMAT(' NEW FRAME')
      END IF
      RETURN
      END
