program filt

c***********************************************************************
c  This is a FORTRAN program to construct square response functions for
c  weights derived from finite Fourier expansions.
c
c  The dimensions of the arrays below have been set so that a maximum
c  number of iwx weights can be used to describe the weighting function
c  This gives a total of 2*iwx-1 weights when you take into account the 
c  assumed symmetry of the weighting function.
c
c  The first input parameter is the position within the Nyquist interval
c  where you want to put the cut.  The Nyquist interval runs from zero
c  to pi in the frequency units used here.  The low-pass cutoff is
c  put at alf*pi, where alf must be on the interval 0<alf<1.
c
c  You can take as many weights (nm) as you want, but the dimension of
c  w(iwx) must be changed if you want more than iwx coefficients.
c
c  The program also lets you decide whether you want to smooth the
c  weights or not.  In general you would always do so.
c************************************************************************

      parameter(ifx=100, iwx=150)
c      ifx = the number of frequencies calculated on the interval 0-pi
c      iwx = the maximum possible number of weights
      character*1 answer
      dimension om(ifx),w(iwx),h(ifx),f(ifx)
      open(12,file='filt.out',form='formatted')
c  First take inputs
      write(*,*) 'where (0-1.0) do you want the cut?'
      read(*,*) alf
      write(*,*) 'alf=',alf
      write(*,*) 'How many weights?'
      read(*,*) nw
      write(*,*) nw,'Weights'
c  Calculate weights
      w(1)=alf
      pi=acos(-1.0)
      api=alf*pi
      do 10 i=2,nw
      k=i-1
      w(i)=sin(api*k)/(pi*k)
   10 continue
c  Smooth weights
c      write(*,*) 'Do you want to Lanzcos filter the weights?'
c      read(*,'(a)') answer
      answer = 'y'
      if(answer.eq.'y') then
      do 45 i=2,nw
      k=i-1.
      xkk=pi*k/float(nw-1)
      w(i)=w(i)*sin(xkk)/xkk
   45 continue
      endif
c      open(22,file='weights',form='formatted')
      write(*,*) 'Here are the weights for alpha=',alf,' nw=',nw
      write(12,*) 'Here are the weights for alpha=',alf,' nw=',nw
      if( answer.eq.'y') then
      write(*,*) 'Lanzcos smoothing has been applied'
      write(12,*) 'Lanzcos smoothing has been applied'
      endif
      write(*,*) '  Weight index,  Weight'
      write(12,*) '  Weight index,  Weight'
      do 11 i=1,nw
      k=i-1
      write(*,*) k,w(i)
      write(12,77) k,w(i)
   77 format(1x,i6,f12.6)
  11  continue
c  Now lets calculate the response function
      xn2=(ifx-1.)*2.0
      pi2=pi*2.0
      do 20 i=1,ifx
      f(i)=(i-1.)/xn2
      om(i)=f(i)*pi2
  20  continue
      do 30 j=1,ifx
      h(j)=w(1)	
      do 25 i=2,nw
      h(j)=h(j)+2.0*w(i)*cos(om(j)*(i-1.))
  25  continue
  30  continue
      write(*,*) 'Here is the response function on zero to pi'
      write(*,*) ' Frequency,   Response'
      write(12,*) 'Here is the response function on zero to pi'
      write(12,*) ' Frequency,   Response'
      iffx=ifx/2
c      do 35 j=1,iffx
      do 35 j=1,ifx
      write(*,78) f(j),h(j)
      write(12,78) f(j),h(j)
   78 format(1x,2f12.6)
   35 continue
      stop
      end