!**********************************************************************!
!
! File: fsaite.f90 (19-May-2000)
!                  (04-May-2001)
!                  (27-May-2005)
!
! Loesung der Wellengleichung fuer die gezupfte Saite
!
!   (d/dt)**2 A(x,t) = c**2 (d/dx)**2 A(x,t)
!
! fuer das Intervall [0,1] und c=1 mit Fourier-Reihe
!
! nx............Anzahl Punkte auf x-Achse
! nk............Anzahl Terme in Fourier-Entwicklung
! t.............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

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

  pi=4.0*atan(1.0)

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

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

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

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

  enddo

  close(unit=iout)

end program fsaite

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