PROGRAM rc
*     simulation of RC circuit with ac voltage source
      IMPLICIT NONE
      REAL*8 I, Q, R, V0, dt, omega, t, tau, tmax
      INTEGER IC1_, IC2_, IR_
      DATA IC1_/1/
      DATA IC2_/2/
      CALL GWopen(IR_, 0)
      CALL initial(Q,R,tau,V0,omega,tmax,t,dt)
      CALL set_up_windows(IC1_,IC2_,REAL(V0),REAL(tmax))
      CALL source_voltage(IC1_,V0,omega,t)
      CALL output_voltage(IC2_,I,R,t)
      DO WHILE(t .LE. tmax)
         CALL scope(Q,I,R,tau,V0,omega,t,dt)
         t = t + dt
         CALL source_voltage(IC1_,V0,omega,t)
         CALL output_voltage(IC2_,I,R,t)
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(Q,R,tau,V0,omega,tmax,t,dt)
      IMPLICIT NONE
      REAL*8 C, f, period, pi, Q, R, tau, V0, omega, tmax, t, dt
      PARAMETER(pi = 3.141592653589793D0)
      t = 0
      V0 = 1                                ! amplitude of external voltage
      WRITE(*,'(A,$)') 'external voltage frequency (Hertz) = '
      READ(*,*) f
      omega = 2*pi*f                        ! angular frequency
      WRITE(*,'(A,$)') 'resistance (ohms) = '
      READ(*,*) R
      WRITE(*,'(A,$)') 'capacitance (farads) = '
      READ(*,*) C
      WRITE(*,'(A,$)') 'time step dt = '
      READ(*,*) dt
      Q = 0
      tau = R*C                             ! relaxation time
      period = 1/f                          ! period of external frequency
      IF(period .GT. tau) THEN
         tmax = 2*period
      ELSE
         tmax = 2*tau
      END IF
      END
C
      SUBROUTINE set_up_windows(IC1_,IC2_,V0,tmax)
      IMPLICIT NONE
      REAL Vmin, tmin, V0, tmax, GWaspect
      INTEGER IR_, IC1_, IC2_
      CALL GWvport(IR_, 0.0, 0.0, GWaspect(1), 0.5)
      tmin = 0
      Vmin = -V0
      CALL draw_axes(tmin,tmax,Vmin,V0)
      CALL GWputtxt(IR_, tmin, V0, 'source voltage')
      CALL GWsavevp(IR_, IC1_)
      CALL GWvport(IR_, 0.0, .5, GWaspect(1), 1.0)
      CALL draw_axes(tmin,tmax,Vmin,V0)
      CALL GWputtxt(IR_, tmin, V0, 'voltage drop across resistor')
      CALL GWsavevp(IR_, IC2_)
      END
C
      SUBROUTINE scope(Q,I,R,tau,V0,omega,t,dt)
*     compute voltage drops
      IMPLICIT NONE
      REAL*8 Q, I, R, tau, V0, omega, t, dt, V
      I = V(V0,omega,t)/R - Q/tau
      Q = Q + I*dt
      END
C
      SUBROUTINE source_voltage(IC1_,V0,omega,t)
      IMPLICIT NONE
      REAL*8 V0, omega, t, V
      REAL x, y
      INTEGER IR_, IC1_
      SAVE x, y
      DATA x, y/2*-100.0/
      CALL GWselvp(IR_, IC1_)
      CALL GWsetpen(IR_, 16, -1, -1, -1)
      CALL GWmove2(IR_, x, y)
      x = REAL(t)
      y = REAL(V(V0,omega,t))
      CALL GWline2(IR_, x, y)
      END
C
      SUBROUTINE output_voltage(IC2_,I,R,t)
      IMPLICIT NONE
      REAL*8 I, R, t
      REAL x, y
      INTEGER IR_, IC2_
      SAVE x, y
      DATA x, y/2*-100.0/
      CALL GWselvp(IR_, IC2_)
      CALL GWsetpen(IR_, 13, -1, -1, -1)
      CALL GWmove2(IR_, x, y)
      x = REAL(t)
      y = REAL(I*R)
      CALL GWline2(IR_, x, y)
      END
C
      FUNCTION V(V0,omega,t)
      IMPLICIT NONE
      REAL*8 V,V0,omega,t
      V = V0*cos(omega*t)
      END