!**********************************************************************!
!
! File: frlkreis.f90 (19-May-2000)
!                    (11-May-2001)
!                    (06-May-2005)
!                    (27-May-2005)
!
! 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

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

  pi=4.0*atan(1.0)

  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
    elseif(modulo(iw,4)==0) then
      a(iw)=-(2.0/(pi*(iw**2-1)))*umax
    else
      a(iw)=0.0
    endif
  enddo

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

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

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

  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)
    enddo

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

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

  enddo

  close(unit=iout)

end program frlkreis

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