PROGRAM DCAPRC
C
C     DISCHARGE OF A CAPACITOR IN AN RC CIRCUIT
C
C     -----------------------------------------
C     AUTHOR: ROBERT EHRLICH (GEORGE MASON U    1973)
C     Revised: Peter Signell (Mich. St. U.      3/86)
C     Revised: Gary Herzenstiel (Mich. St. U.   3/86)
C     Revised: E.J.D.Kales (Mich. St. U.       11/87)
C
C     Calculate and plot the amount of charge remaining on a
C     discharging capacitor as a function of time.
C     Calculate an approximate value for the charge using an
C     iterative procedure in which the change in the charge
C     during a time interval DT is given by DQ= -Q*DT/TAU.
C     Also caluclate exact value of the charge using
C     function ECP(-1/TAU).
C
C     SUBROUTINES:
C          GETINF  - GETS VALUES FOR TAU,DT,N FROM USER
C          CALC    - CALCULATES CHARGE REMAINING AND
C                    PUTS IT INTO GRAPHICAL FORMAT
C          DISPLA  - DISPLAYS THE GRAPH
C          QUERY   - ASKS USER IF TO RUN THE PROGRAM AGAIN
C
C                       D C A P R C C
C                             :
C         -----------------------------------------------
C         :            :               :                :
C      GETINF        CALC            DISPLA           QUERY
C     (TAU,DT,N1)  (TAU,DT,N1,      (N1,X,TARRAY,    (AGAIN)
C                  TARRAY,Q1ARRY,   Q1ARRY,
C                  Q2ARRY,X)       Q2ARRY,X)
C
C     ***** DECLARATIONS *****
      REAL N,TAU,DT
      REAL TARRAY(501),Q1ARRY(501),Q2ARRY(501)
      INTEGER N1
      CHARACTER X(501)*101,AGAIN
      OPEN(UNIT=3,FILE='m350p1f.out')
C                    
C     ***** M A I N   B O D Y *****
C
 10   WRITE(3,'(A1//)') 12
      CALL GETINF (                            TAU,DT,N1)
      CALL CALC   (TAU,DT,N1,   TARRAY,Q1ARRY,Q2ARRY,X)
      CALL DISPLA (N1,X,TARRAY,Q1ARRY,Q2ARRY          )  
      CALL QUERY  (                                AGAIN)
      IF(AGAIN.EQ.'Y') GOTO 10
C     IF USER IS DONE, STOP
      STOP
      END
C     ***** M A I N   P R O G R A M    E N D *****
C
C     ***** G E T I N F *****
C
      SUBROUTINE GETINF(   TAU,DT,N1)
C
C     GETS VALUES FOR TAU,DT,N.  N1 IS INTEGER FORM OF N. 
C
C     ***** DECLARATIONS *****
      REAL TAU,DT,N
      INTEGER N1
C
C     TAU = time constant
C     DT  = time step
C     N   = number of time steps
C
      PRINT *
      PRINT *,'Please supply TAU, DT & N values:'
      WRITE(*,20) 'TAU='
      READ *,TAU
      WRITE(*,20) ' DT='
      READ *,DT
      WRITE(*,20) '  N='
      READ *,N
  20  FORMAT(A4/)
  21  FORMAT(F8.5)
      N1=N
      PRINT *
      WRITE(*,34) TAU,DT,N1
      WRITE(3,34) TAU,DT,N1
 34   FORMAT('m350p1f.for: Discharge of a Capacitor:'//
     +' TAU=',F8.5,'  DT=',F8.5,'   N=',I4/)
      RETURN                                         
      END
C
C     ***** C A L C *****
C
      SUBROUTINE CALC(TAU,DT,N1,   TARRAY,Q1ARRY,Q2ARRY,X)
C
C     find amount of charge remaining and put into X array 
C
C     ***** DECLARATIONS *****
      REAL TAU,DTM1,DQ,Q1,Q2,T
      REAL TARRAY(501),Q1ARRY(501),Q2ARRY(501)
      INTEGER J,M1,M2,N1
      CHARACTER X(501)*101
C
C     Q1 = exact value of the charge remaining
C     Q2 = approximate value of the charge remaining
C     T  = time
C     X  = graph array for Display
C
C     at time T=0, both Q1 and Q2 are assumed to be 1.0
C
      T=0.0
      Q1=1.0
      Q2=1.0
C     loop over rows
      DO 30,J=1,N1
C
C          Initialize & set boundaries for graph:
           X(J)='.                                                      
     +    .'
C
C          compute M1, between 1 and 80 (or 101)
C          corresponding to vaules of Q1, which can range
C          between 0.0 and 1.0, and likewise M2 for Q2.
C
CC large or small:
           M1=59.0*Q1+1.0
           M2=59.0*Q2+1.0
CC           M1=100.0*Q1+1.0
CC           M2=100.0*Q2+1.0
C
C          store '*' to represent Q1, '+' to represent Q2
C
           X(J)(M1:M1)='*'
           X(J)(M2:M2)='+'
C
C          store values of T,Q1,Q2 in arrays for later
C          printing
C
           TARRAY(J)=T
           Q1ARRY(J)=Q1
           Q2ARRY(J)=Q2
C
C          calculate the next values of Q1 and Q2
C
           T=T+DT
           DQ=-Q2*DT/TAU
           Q2=Q2+DQ
           Q1=EXP(-T/TAU)
 30   CONTINUE
      RETURN
      END
C
C     ***** D I S P L A *****
C
      SUBROUTINE DISPLA(N1,X,TARRAY,Q1ARRY,Q2ARRY     )
C
C     ***** DECLARATIONS *****
      REAL TARRAY(501),Q1ARRY(501),Q2ARRY(501)
      INTEGER J,N1
      CHARACTER X(501)*101
C
      WRITE(3,22)
22    FORMAT('  TIME  EXACT APPROX'/
     +       '         Q     Q   0'/
     +       '         *     +    ')
      WRITE(3,23)
23    FORMAT(19X,41H.........................................,
CC large or small:
     +19H...................)
CC   +60H...................................................
CC   +.........)
C
C     loop to display graph
C
CC  for large printer:
CC      WRITE(3,24) (TARRAY(J),Q1ARRY(J),
CC   +      Q2ARRY(J),X(J),J=1,N1)
CC 24   FORMAT(1X,F5.2,2F6.3,1X,A101)
CC  for small printer:
      WRITE(3,24) (TARRAY(J),Q1ARRY(J),Q2ARRY(J),X(J)(1:60),J=1,N1)
 24   FORMAT(1X,F5.2,2F6.3,1X,A60)
      WRITE(3,23)
      RETURN
      END
C
C     ***** Q U E R Y *****
C
      SUBROUTINE QUERY(    AGAIN)
C 
C     ***** DECLARATIONS *****
      CHARACTER AGAIN
C
      PRINT *
      WRITE(*,62)
 62   FORMAT(
     +'Would you like to do this again using new values? (Y/N)? ')
      READ( *,'(A1)')AGAIN
      RETURN
      END