!**********************************************************************!
!
! File: rlckreis-rk4.f90 (24-May-2001)
!                        (03-Jun-2002)
!                        (27-May-2005)
!
! Einschwingverhalten eines RLC-Kreises
!
!   dI/dt = (U - R*I - Q/C)/L
!
!   dQ/dt = I
!
! mit Anfangsbedingung I(t_0)=0, Q(t_0)=0 und
!
!   U(t) = U_max*cos(w*t)   fuer cos(w*t) > 0
!
!          0                sonst
!
! Integration mit Runge-Kutta Verfahren 4. Ordnung
!
!   R....................Widerstand
!   L....................Induktivitaet
!   C....................Kapazitaet
!   U_max................Spannung (Maximalamplitude)
!   w....................Kreisfrequenz
!   dt...................Zeitschritt
!   t_0,t_max............Anfangs/Endzeitpunkt
!   "rlckreis-rk4.dat"...Ausgabefile
!
!**********************************************************************!

module rlcrk4_m

  implicit none

  public::f,rk4
  real,public::c,r,umax,w,xl

contains

!**********************************************************************!

  subroutine f(dxidt,dqdt,xi,q,t)

!**********************************************************************!

    real,intent(in)::q,t,xi
    real,intent(out)::dqdt,dxidt
    real::u

    if(cos(w*t)>0.0) then
      u=umax*cos(w*t)
    else
      u=0.0
    end if

    dxidt=(u-r*xi-q/c)/xl
    dqdt=xi

    return

  end subroutine f

!**********************************************************************!

  subroutine rk4(u1n,u2n,tn,u1,u2,t)

!**********************************************************************!

    real,intent(in)::t,tn,u1,u2
    real,intent(out)::u1n,u2n
    real::dt,f11,f12,f21,f22,f31,f32,f41,f42

    dt=tn-t

    call f(f11,f12,u1,u2,t)
    call f(f21,f22,u1+0.5*dt*f11,u2+0.5*dt*f12,t+0.5*dt)
    call f(f31,f32,u1+0.5*dt*f21,u2+0.5*dt*f22,t+0.5*dt)
    call f(f41,f42,u1+dt*f31,u2+dt*f32,t+dt)

    u1n=u1+dt*(f11+2.0*f21+2.0*f31+f41)/6.0
    u2n=u2+dt*(f12+2.0*f22+2.0*f32+f42)/6.0

    return

  end subroutine rk4

end module rlcrk4_m

!**********************************************************************!

program rlckreis_rk4

!**********************************************************************!

  use rlcrk4_m
  
  implicit none

  integer,parameter::iout=2

  integer::it
  real::dt,q,qn,t,t0,tmax,tn,xi,xin

  write(unit=*,fmt="(""          R(Ohm)="")",advance="no")
  read(unit=*,fmt=*) r
  write(unit=*,fmt="(""           L(Hy)="")",advance="no")
  read(unit=*,fmt=*) xl
  write(unit=*,fmt="(""            C(F)="")",advance="no")
  read(unit=*,fmt=*) c
  write(unit=*,fmt="(""        U_max(V)="")",advance="no")
  read(unit=*,fmt=*) umax
  write(unit=*,fmt="(""     w/2*pi(1/s)="")",advance="no")
  read(unit=*,fmt=*) w
  w=8.0*atan(1.0)*w
  write(unit=*,fmt="("" t_0,dt,t_max(s)="")",advance="no")
  read(unit=*,fmt=*) t0,dt,tmax

  open(unit=iout,file="rlckreis-rk4.dat",status="replace", &
       action="write",form="formatted")

  write(unit=iout,fmt="("" #"")")
  write(unit=iout,fmt="("" #   t(s)        I(A)           Q(C)"")")
  write(unit=iout,fmt="("" #"")")

  t=t0
  xi=0.0
  q=0.0
  write(unit=iout,fmt="(tr1,f10.5,2(tr1,es14.7))") t,xi,q

  do it=1,int((tmax-t0)/dt+0.5)

    tn=t0+it*dt
    call rk4(xin,qn,tn,xi,q,t)

    t=tn
    xi=xin
    q=qn
    write(unit=iout,fmt="(tr1,f10.5,2(tr1,es14.7))") t,xi,q

  end do

  close(unit=iout)

end program rlckreis_rk4

!**********************************************************************!