      program sp511
      implicit none
      integer j
      real*8 f, Imax

      f = 2d2
      call sp510( _, ____ )		  ! Hz,  Maximum current
      write(*,*) 'f =', f, '  Imax =', Imax

      open(1,file='sp510.dat')
      do j=10,60
        f = 10**(j/10.)			  ! Hz
        call sp510( _, ____ )		  ! Hz,  Maximum current
        write(*,'(1p2e12.4)') f, Imax
        write(1,'(1p2e12.4)') f, Imax
      enddo
      close(1)
      stop
      end

      subroutine sp510( _, ____ )	  ! Hz,  Maximum current
      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 )    ! volt
      data t,Q,I/0,0,0/     ! s:second, C:coulomb, A:ampere
      Imax  = 0
      if( 1d0.le.f .and. f.lt.3d0 ) Nsample = 30000
      if( 3d0.le.f .and. f.lt.5d0 ) Nsample = 10000
      if( 5d0.le.f .and. f.lt.1d1 ) Nsample =  5000
      if( 1d1.le.f .and. f.lt.3d1 ) Nsample =  3000
      if( 3d1.le.f .and. f.lt.1d2 ) Nsample =  1000
      if( 1d2.le.f .and. f.lt.1d3 ) Nsample =   500
      if( 1d3.le.f                ) Nsample =   200
      period= 1/f
      dt    = period / Nsample
      jmax  = period/dt*6
      jmax2 = period/dt*2
      do j=1,jmax
        Vs    = ext_volt_cos      (t,f,V0)           ! sp401.f
        A     = - (R/L)*I - Q/(L*C) + Vs/L           ! 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
        Amid  = - (R/L)*Imid - Qmid/(L*C) + Vsmid/L  ! A=dI/dt
        I     = I + Amid*dt
        Q     = Q + Imid*dt
        VR    = I*R
        VC    = Q/C
        VL    = Amid*L
        VAL   = Vs - (VR+VC)
        t     = t + dt
!       write( *, '(1p8e12.4)' ) t,Vs,I,VR,VC,VL,VAL,Q
        if( j.gt.jmax2 .and. I.gt.Imax ) Imax = I
      enddo
      return
      end
