!**********************************************************************!
!
! File: rlkreis.f90 (26-May-2008)
!
! Einschwingverhalten eines RL-Kreises
!
!   L*dI/dt + R*I = U,     I(t_0)=0
!
! bei Anregung mit
!
!          U_max*cos(w*t)   fuer cos(w*t) > 0
!   U(t) =
!          0                sonst
!
! Integration mit Verfahren 2. Ordnung mit Hilfsschritt
! (Runge-Kutta Verfahren 2. Ordnung)
!
! R...............Widerstand
! L...............Induktivitaet
! U_max...........Spannung (Amplitude)
! w...............Kreisfrequenz
! dt..............Zeitschritt
! t_0.............Startzeitpunkt
! t_max...........Endzeitpunkt
!
! "rlkreis.dat"...Ausgabefile
!
!**********************************************************************!

module rlkreis_m

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

  implicit none

  private

  integer,parameter,public::iout=2

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

  public::f

contains

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

  function f(xi,t) result(y)

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

    real,intent(in)::t,xi
    real::y

    if(cos(w*t)>=0.0) then
      y=(-r*xi+umax*cos(w*t))/xl
    else
      y=(-r*xi)/xl
    end if

    return

  end function f

end module rlkreis_m

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

program rlkreis

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

  use rlkreis_m

  implicit none

  integer::it,itmax
  real::dt,pi,t,t0,t1,tmax,tn,xi,xi1,xin

  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") "        U_max[V]="
  read(unit=*,fmt=*) umax
  write(unit=*,fmt="(a)",advance="no") "   w/(2*pi)[1/s]="
  read(unit=*,fmt=*) w
  write(unit=*,fmt="(a)",advance="no") " t_0,dt,t_max[s]="
  read(unit=*,fmt=*) t0,dt,tmax
  w=2.0*pi*w

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

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

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

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

  do it=1,itmax

    t1=t+0.5*dt           !Hilfsschritt
    xi1=xi+0.5*dt*f(xi,t)

    tn=t0+it*dt           !Vollschritt
    xin=xi+dt*f(xi1,t1)

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

  end do

  close(unit=iout)

end program rlkreis

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