!**********************************************************************!
!
! File: heavi.f90 (19-May-2000)
!                 (24-Apr-2006)
!                 (11-May-2007)
!                 (25-Mar-2008)
!
! Fouriersynthese der Heaviside-Funktion
!
!          0 fuer x < 0
!   H(x) =
!          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

  integer,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

  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

  dx=(xmax-xmin)/nx
  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)
  end do

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

  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)
    end do

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

  end do

  close(unit=iout)

end program heavi

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