!**********************************************************************!
!
! File: fsaite.f90 (19-May-2000)
!                  (28-Apr-2006)
!                  (18-May-2007)
!                  (23-Mar-2008)
!
! Loesung der Wellengleichung fuer die gezupfte Saite
!
!   (d/dt)**2 A(x,t) = c**2 (d/dx)**2 A(x,t)
!
! mit c=1 und Maximalauslenkung A=1 auf dem Intervall [0,1]
! mittels Fourier-Reihe
!
! nx.............Anzahl Punkte auf x-Achse
! nk.............Anzahl Terme in Fourier-Entwicklung
! t..............gewuenschter Zeitpunkt
!
! "fsaite.dat"...Ausgabefile
!
!**********************************************************************!

program fsaite

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

  implicit none

  integer,parameter::iout=2

  integer::ik,ix,nk,nx
  real::pi,s,t,x
  real,dimension(:),allocatable::b

  pi=4.0*atan(1.0)

  write(unit=*,fmt="(a)",advance="no") " nx="
  read(unit=*,fmt=*) nx
  write(unit=*,fmt="(a)",advance="no") " nk="
  read(unit=*,fmt=*) nk
  write(unit=*,fmt="(a)",advance="no") "  t="
  read(unit=*,fmt=*) t

  allocate(b(1:nk))

  do ik=1,nk,2
    if(modulo(ik,4)==1) then
      b(ik)=8.0/(ik*pi)**2
    else
      b(ik)=-8.0/(ik*pi)**2
    end if
  end do

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

  do ix=0,nx

    x=real(ix)/nx
    s=0.0
    do ik=1,nk,2
      s=s+b(ik)*cos(ik*pi*t)*sin(ik*pi*x)
    end do

    write(unit=iout,fmt="(tr1,f8.5,tr1,es14.7)") x,s

  end do

  close(unit=iout)

end program fsaite

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