!**********************************************************************!
!
! File: frlkreis.f90 (19-May-2000)
!                    (05-May-2006)
!                    (24-Mar-2008)
!
! Strom in einem RL-Kreis
!
!   L*dI/dt + R*I = U
!
! bei periodischer Anregung mit positiven Cosinus-Halbwellen
! (Periode 1(s))
!
! Loesung mit Fourier-Ansatz ueber dem Intervall [-0.5(s),0.5(s)]
!
! nt...............Anzahl Punkte auf t-Achse
! nw...............Anzahl Terme in Fourier-Entwicklung (w=omega)
! R................Widerstand
! L................Induktivitaet
! U_max............Spannung (Amplitude)
!
! "frlkreis.dat"...Ausgabefile
!
!**********************************************************************!

program frlkreis

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

  implicit none

  integer,parameter::iout=2

  integer::it,iw,nt,nw
  real::pi,r,su,t,umax,xl
  real,dimension(:),allocatable::a
  complex::si,z
  complex,dimension(:),allocatable::c

  pi=4.0*atan(1.0)

  write(unit=*,fmt="(a)",advance="no") "       nt="
  read(unit=*,fmt=*) nt
  write(unit=*,fmt="(a)",advance="no") "       nw="
  read(unit=*,fmt=*) nw
  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

  allocate(a(0:nw),c(-nw:nw))

  a(0)=(2.0/pi)*umax
  a(1)=0.5*umax
  do iw=2,nw
    if(modulo(iw,4)==2) then
      a(iw)=(2.0/(pi*(iw**2-1)))*umax
    else if(modulo(iw,4)==0) then
      a(iw)=-(2.0/(pi*(iw**2-1)))*umax
    else
      a(iw)=0.0
    end if
  end do

  do iw=0,nw
    z=cmplx(r,2.0*pi*iw*xl)
    c(iw)=0.5*a(iw)/z
    c(-iw)=conjg(c(iw))
  end do

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

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

  do it=0,nt

    t=-0.5+real(it)/nt

    su=0.5*a(0)
    do iw=1,nw
      su=su+a(iw)*cos(2.0*pi*iw*t)
    end do

    si=0.0
    do iw=-nw,nw
      si=si+c(iw)*exp(cmplx(0.0,2.0*pi*iw*t))
    end do

    write(unit=iout,fmt="(tr1,f8.5,2(tr1,es14.7))") t,su,real(si)

  end do

  close(unit=iout)

end program frlkreis

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