!**********************************************************************!
!
! File: heavi.f90 (19-May-2000)
!                 (04-May-2001)
!                 (27-May-2005)
!
! Fouriersynthese der Heaviside-Funktion
!
!   H(x) = 0 fuer x < 0
!          1 fuer x > 0
!
! auf dem Intervall [-0.5,0.5]
!
! nx............Anzahl Punkte auf x-Achse
! nk............Anzahl Terme in Fourier-Entwicklung
! "heavi.dat"...Ausgabefile
!
!**********************************************************************!

program heavi

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

  implicit none

  real,parameter::iout=2
  real,parameter::xmax=0.5,xmin=-xmax

  integer::ik,ix,nk,nx
  real::a0,dk,dx,pi,s,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

  dx=(xmax-xmin)/nx
  pi=4.0*atan(1.0)
  dk=2.0*pi/(xmax-xmin)

  allocate(b(1:nk))

  a0=1.0
  do ik=1,nk,2
    b(ik)=(4.0/(ik*dk))/(xmax-xmin)
  enddo

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

  do ix=0,nx

    x=xmin+ix*dx
    s=0.5*a0
    do ik=1,nk,2
      s=s+b(ik)*sin(ik*dk*x)
    enddo

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

  enddo

  close(unit=iout)

end program heavi

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