      PROGRAM RLCCIR
C
C     RLC SERIES 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 current as a function of time for an RLC 
C     series circuit in which the capacitor is initially charged,
C     and the current is initially zero.
C
C     SUBROUTINES:
C          GETINF  - get values for R,L,C,DT,N,V0,W
C          CALC    - calculate current and put in graph
C          DISPLA  - display the graph 
C          QUERY   - ask if user wants to repeat
C
C                            R L C S C I R
C                                  :
C            ----------------------------------------------
C            :            :                  :            :
C         GETINF        CALC              DISPLA        QUERY
C       (R,L,C,DT,  (R,L,C,DT,N1,V0,    (N1,X,TARRAY,  (AGAIN)
C        N1,V0,W)   W,X,TARRAY,IARRAY)    IARRAY)
C                                           
C     ***** DECLARATIONS *****
      REAL R,L,C,DT,V0,W,TARRAY(501),IARRAY(501)
      INTEGER N1
      CHARACTER X(501)*101,AGAIN
      OPEN(UNIT=3,FILE='(T:83)m351p1f.out')
CC      OPEN(UNIT=3,FILE='(T:120)OUTPUT.LST')
C                    
C     ***** M A I N   B O D Y *****
C
 10   WRITE(3,'(A1)') 12
      CALL GETINF (                   R,L,C,DT,N1,V0,W)
      CALL CALC   (R,L,C,DT,N1,V0,W,   TARRAY,IARRAY,X)
      CALL DISPLA (N1,X,TARRAY,IARRAY                 )  
      CALL QUERY  (                              AGAIN)
      IF(AGAIN.EQ.'Y') GOTO 10
      STOP
      END
C
C     ***** G E T I N F *****
C
      SUBROUTINE GETINF(   R,L,C,DT,N1,V0,W)
C 
C     get values for R,L,C,DT,N,V0, & W.  N1 is integer of N. 
C
C     ***** DECLARATIONS *****
      REAL R,L,C,DT,N,V0,W
      INTEGER N1
C
C     R   = resistance
C     L   = inductance
C     C   = capacitance
C     DT  = time step
C     N   = number of time steps
C     V0  = amplitude of applied voltage
C     W   = frequency of applied voltage
C
      PRINT *,'Please enter R, L, C, DT, N, V0 & W values:'
      WRITE(*,20) ' R='
      READ *,R
      WRITE(*,20) ' L='
      READ *,L
      WRITE(*,20) ' C='
      READ *,C
      WRITE(*,20) 'DT='
      READ *,DT 
      WRITE(*,20) ' N='
      READ *,N 
      WRITE(*,20) 'V0='
      READ *,V0
      WRITE(*,20) ' W='
      READ *,W
  20  FORMAT(A4/)
      N1=N
      WRITE(*,34) R,L,C,DT,N1,V0,W
      WRITE(3,34) R,L,C,DT,N1,V0,W
 34   FORMAT('M351P1.FOR'/
     + 'RLC Series Circuit:    R=',F8.4/20X,'   L=',F8.4/
     +20X,'   C=',F8.4/20X,'  DT=',F8.4/20X,'   N=',I4/
     +20X,'  V0=',F8.4/20X,'   W=',F8.4/)
      RETURN                                         
      END
C
C     ***** C A L C *****
C
      SUBROUTINE CALC(R,L,C,DT,N1,V0,W,   TARRAY,IARRAY,X)
C 
C     find the current for an RLC using N1 intervals
C
C     ***** DECLARATIONS *****
      REAL R,L,C,DT,V0,W,T,Q,I,I1,VA,DIDT,VA1,DIDT1
      REAL TARRAY(501),IARRAY(501)
      INTEGER J,M,N1
      CHARACTER X(501)*101
C
C     T  = time
C     X  = graph array
C
C     at time T=0, I=0 and Q=1
C
      T=0.0
      Q=1.0
      I=0.0
C
      DO 30,J=1,N1
C
C          initialize graph line and set boundaries:
CC small or large:
           X(J)='.                               '
     +       //'|                               .'
CC           X(J)='.                                                 '
CC     +       //'|                                                 .'
C
C          compute an index M, in the range 0 to 65 or 101 to correspond
C          to a value of the current in the range -1.0 to +1.0
C
CC small or large:
           M=32*(I+1.03)
CC           M=50.0*(I+1.03)
           IF(M.LE.0) M=1
CC small or large:
           IF(M.GT.80) M=65
CC           IF(M.GT.100) M=101
C
C          store '*' in X to represent M
C
           X(J)(M:M)='*'
C
C          store values of T & I into arrays for later printing
C
           TARRAY(J)=T
           IARRAY(J)=I
C
C          calculate the next values 
C
           VA=V0*SIN(W*T)
C
C          DIDT is the time interval of I
C
           DIDT=(1.0/L)*(VA-I*R-Q/C)
C
C          approximate new Q and I from old values
C
           Q1=Q+I*DT
           I1=I+DIDT*DT
           VA1=V0*SIN(W*(I+DT))
           DIDT1=(1.0/L)*(VA1-I1*R-Q1/C)
C
C          get improved values for Q and I
C
           Q=Q+(I+I1)*DT/2.0
           I=I+(DIDT+DIDT1)*DT/2.0
           T=T+DT
30    CONTINUE
      RETURN
      END
C
C     ***** D I S P L A *****
C
      SUBROUTINE DISPLA(N1,X,TARRAY,IARRAY    )
C 
C     ***** DECLARATIONS *****
      REAL TARRAY(501),IARRAY(501)
      INTEGER J,N1
      CHARACTER X(501)*101
C
      WRITE(3,52)
 52   FORMAT('  TIME    CURRENT')
      WRITE(3,53)
 53   FORMAT(18X,41H.........................................,
CC for large printer:
CC     +60H............................................................)
CC      WRITE(3,54) (TARRAY(J),IARRAY(J),X(J),J=1,N1)
CC 54    FORMAT(F7.2,F9.4,2X,A101)
CC for small printer:
     +24H........................)
      WRITE(3,54) (TARRAY(J),IARRAY(J),X(J)(1:65),J=1,N1)
 54   FORMAT(F7.2,F9.4,2X,A65)
      WRITE(3,53)
      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 points? (Y/N)? ')
      READ( *,'(A1)')AGAIN
      RETURN
      END
C     ***** P R O G R A M    E N D *****