!**********************************************************
!*                    PROGRAM CIRCUIT                     *
!* Use Runge-Kutta method to solve a LRC circuit or equi- *
!* valent damped mass-spring problem.                     *
!* ------------------------------------------------------ *
!* SAMPLE RUN:                                            *
!*                                                        *
!* (Solve mass-spring system with:                        *
!*  F=75 N, M=50 kg, K=100 N (C=0.01), D=0.05)            *
!*                                                        *
!* Specify [o]scillating or [n]o oscillating term...      *
!* o                                                      *
!* Give V or F, L or M, R or D, C or 1/K (4L/C-R^2 > 0):  *
!* 75 50 0.05 0.01                                        *
!* one period at t=   4.442879    s                       *
!* Give t_final and number of points:                     *
!* 10 10                                                  *
!*                                                        *
!*     time        q or d    i or speed  analytic q or d  *
!* 1.000000E+00  7.5037E-01  1.4993E+00    6.3310E-01     *
!* 2.000000E+00  2.2485E+00 -1.7480E-03    1.4628E+00     *
!* 3.000000E+00  7.4825E-01 -2.9955E+00    1.0890E+00     *
!* 4.000000E+00 -2.2440E+00  5.4882E-03    1.4358E-01     *
!* 5.000000E+00  7.5549E-01  5.9850E+00    2.2231E-01     *
!* 6.000000E+00  6.7320E+00 -1.4953E-02    1.1913E+00     *
!* 7.000000E+00  7.3505E-01 -1.1958E+01    1.4147E+00     *
!* 8.000000E+00 -1.1202E+01  3.7843E-02    5.1634E-01     *
!* 9.000000E+00  7.8782E-01  2.3892E+01    1.3089E-02     *
!* 1.000000E+01  2.4630E+01 -9.1529E-02    7.5371E-01     *
!* Press any key to continue                              *
!* (For a better accuracy, increase number of points).    *
!* ------------------------------------------------------ *
!* Reference: "Problem Solving with Fortran 90 By David   *
!*             R. Brooks, Springer-Verlag New York, 1997" *
!********************************************************** 
! 
      MODULE RungeKutta
!
      CONTAINS
!------------------------------------
      REAL FUNCTION AofT(x,v,D,K,M,F)
! Calculate acceleration for "mass and spring" problem.
      IMPLICIT NONE
      REAL, INTENT(IN) :: x,v,D,K,M,F
      AofT=-(-F+D*v+K*x)/M
      END FUNCTION AofT
!--------------------------
      SUBROUTINE MassAndSpring(x,v,D,K,M,F,dt)
! Calculate motion for mass and spring problem with constant force term.
      IMPLICIT NONE
      REAL,INTENT(INOUT) :: x,v
      REAL, INTENT(IN) :: D,K,M,F,dt
      REAL k1_x,k2_x,k3_x,k4_x,k1_v,k2_v,k3_v,k4_v
!
! Runge-Kutta coefficients...
      k1_x=v
      k2_x=v+AofT(x,v,D,K,M,F)*dt/2.
      k3_x=v+AofT(x,v,D,K,M,F)*dt/2.
      k4_x=v+AofT(x,v,D,K,M,F)*dt
      k1_v=AofT(x,v,D,K,M,F)
      k2_v=AofT(x+v*dt/2.,v+k1_v*dt/2.,D,K,M,F)
      k3_v=AofT(x+v*dt/2.,v+k2_v*dt/2.,D,K,M,F)
      k4_v=AofT(x+v*dt,v+k3_v*dt,D,K,M,F)
! Propagate solution...
      x=x+(k1_x+2.*k2_x+2.*k3_x+k4_x)*dt/6.
      v=v+(k1_v+2.*k2_v+2.*k3_v+k4_v)*dt/6.
      END SUBROUTINE MassAndSpring
!---------------------------------
      END MODULE RungeKutta
!==========================
      PROGRAM Circuit
!
! Use Runge-Kutta method to solve LRC circuit problems.
! Variable equivalences with mass-and-spring problem:
! V   => force F
! q   => displacement x
! i   => velocity v
! L   => mass m
! R   => damping constant D
! 1/C => spring constant K
      USE RungeKutta, ONLY : MassAndSpring
      IMPLICIT NONE
      REAL i,q,L,C,R,V,t,dt,t_final
      REAL k1_i,k2_i,k3_i,k4_i
      INTEGER j,n
      CHARACTER*1 choice
!
! Choose circuit type...
!
      PRINT *,' Specify [o]scillating or [n]o oscillating term...'
      READ *,choice
      SELECT CASE (choice)
      CASE ('o','O')
! Ld^2q/dt^2+Rdq/dt+q/C=V
10      PRINT *,' Give V or F, L or M, R or D, C or 1/K (4L/C-R^2 > 0):'
          READ *,V,L,R,C
          PRINT *,' one period at t=',4*3.14159*L/SQRT(4.*L/C-R*R),' s'
          PRINT *,' Give t_final and number of points:'
          READ *,t_final,n
        IF ((4.*L/C-R*R) <= 0.) GO TO 10
        q=0.; i=C*V*R/2./L; t=0. !Initial values
        dt=t_final/REAL(n)
        PRINT *,'    time        q or d    i or speed  analytic q or d'
        DO j=1,n
          CALL MassAndSpring(q,i,R,1./C,L,V,dt)
          t=t+dt
          PRINT 1000,t,q,i,C*V*(1.-EXP(-R*t/2./L)*COS(SQRT(4.*L/C-R*R)*t/2./L))
        END DO
1000    FORMAT(1x,es12.6,2es12.4,es14.4)
!
      CASE ('n','N')
! Ldi/dt+Ri=V (no oscillating term)...
        PRINT *,' Give V, L, R:'
        READ *,V,L,R
        PRINT *,' time constant at t=',L/R,' s'
        PRINT *,' Give t_final and number of points:'
        READ *,t_final,n
        i=0.; t=0. !Initial values
        dt=t_final/REAL(n)
        PRINT *,'       time          i        analytic 1'
        DO j=1,n
! Runge-Kutta coefficients...
          k1_i=(V-R*i             )/L
          k2_i=(V-R*(i+k1_i*dt/2.))/L
          k3_i=(V-R*(i+k2_i*dt/2.))/L
          k4_i=(V-R*(i+k3_i*dt)   )/L
! Propagate solution...
          i=i+(k1_i+2.*k2_i+2.*k3_i+k4_i)*dt/6.
          t=t+dt
          PRINT 1000,t,i,V/R*(1.-EXP(-R*t/L))
        END DO
!
      CASE DEFAULT
        PRINT *,' No such choice.  Try again...'
      END SELECT
!
      END

!end of file circuit.f90