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