!**********************************************************************!
!
! File: rlckreis.f90 (24-May-2001)
!                    (18-May-2006)
!                    (27-Mar-2008)
!
! 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 Spannungsquelle
!
!          U_max*cos(w*t)   fuer cos(w*t) > 0
!   U(t) =
!          0                sonst
!
! Integration mit Runge-Kutta Verfahren 4. Ordnung
!
!   R................Widerstand
!   L................Induktivitaet
!   C................Kapazitaet
!   U_max............Spannung (Amplitude)
!   w................Kreisfrequenz
!   dt...............Zeitschritt
!   t_0,t_max........Anfangs/Endzeitpunkt
!
!   "rlckreis.dat"...Ausgabefile
!
!**********************************************************************!

module rlckreis_m

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

  implicit none

  private

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

  public::f,rk4

contains

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

  subroutine f(dudt,u,t)

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

    real,intent(in)::t
    real,dimension(0:),intent(in)::u
    real,dimension(0:),intent(out)::dudt

    real::dqdt,dxidt,q,v,xi

    xi=u(0)
    q=u(1)

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

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

    dudt(0)=dxidt
    dudt(1)=dqdt

    return

  end subroutine f

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

  subroutine rk4(un,u,t,dt)

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

    real,intent(in)::dt,t
    real,dimension(0:),intent(in)::u
    real,dimension(0:),intent(out)::un

    integer::i,n
    real,dimension(0:size(u)-1)::f1,f2,f3,f4

    n=size(u)

    call f(f1,u,t)

    do i=0,n-1
      un(i)=u(i)+0.5*dt*f1(i)
    end do
    call f(f2,un,t+0.5*dt)

    do i=0,n-1
      un(i)=u(i)+0.5*dt*f2(i)
    end do
    call f(f3,un,t+0.5*dt)

    do i=0,n-1
      un(i)=u(i)+dt*f3(i)
    end do
    call f(f4,un,t+dt)

    do i=0,n-1
      un(i)=u(i)+dt*(f1(i)+2.0*(f2(i)+f3(i))+f4(i))/6.0
    end do

    return

  end subroutine rk4

end module rlckreis_m

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

program rlckreis

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

  use rlckreis_m
  
  implicit none

  integer,parameter::iout=2

  integer::it,itmax
  real::dt,pi,q,t,t0,tmax,xi
  real,dimension(0:1)::u,un

  pi=4.0*atan(1.0)

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

  itmax=int((tmax-t0)/dt+0.5)

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

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

  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,itmax

    u(0)=xi
    u(1)=q

    call rk4(un,u,t,dt)

    t=t0+it*dt
    xi=un(0)
    q=un(1)

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

  end do

  close(unit=iout)

end program rlckreis

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