      program sp500
      implicit none
      real*8 pi
      parameter ( pi=3.14159 )      
      real*8 Imax                         ! Maximum current
      real*8 R, C, L
      parameter( R=100, C=3d-6, L=2d-3 )  ! ohm, farad:F, henry:H
      integer j, jmax, jmax2, Nsample
      real*8 t,dt, Q,I,A, Qmid,Imid,Amid, VR,VC,VL,VAL, tmid,Vsmid
      real*8 V0, f, period, Vs, ext_volt_cos, ext_volt_rectangle
      parameter( V0=1, f=2d2 )    ! volt, Hz
      data t,Q,I,Imax/0,0,0,0/     ! s:second, C:coulomb, A:ampere
      Nsample= 1000
      period = 1/f
      dt     = period / Nsample
      write(*,*)      'f =', f,      '  Nsample =', Nsample
      write(*,*) 'period =', period,      '  dt =', dt
      jmax   = period/dt*6
      jmax2  = period/dt*2
      open( 1, file='sp500.dat' )
      do j=1,jmax
        Vs    = ext_volt_cos      (t,f,V0)           ! sp401.f
!       Vs    = ext_volt_rectangle(t,f,V0)           ! sp402.f
        A     = - _______ - _______ + ____           ! A=dI/dt
        Imid  = I + A*dt/2
        Qmid  = Q + I*dt/2
        tmid  = t + dt/2
        Vsmid = ext_volt_cos      (tmid,f,V0)        ! sp401.f
!       Vsmid = ext_volt_rectangle(tmid,f,V0)        ! sp402.f
        Amid  = - (R/L)*Imid - Qmid/(L*C) + Vsmid/L  ! A=dI/dt
        I     = I + ____*dt
        Q     = Q + ____*dt
        VR    = I*R
        VC    = Q/C
        VL    = Amid*L
        VAL   = Vs - (VR+VC)
        t     = t + dt
        write( 1, '(1p8e12.4)' ) t,Vs,I,VR,VC,VL,VAL,Q
        if( j.gt.jmax2 .and. I.gt.Imax ) Imax = I
      enddo
      write(*,*) 'Imax =', Imax
      close(1)
      stop
      end
