c
c Illustrations for Lanczos Filter as discussed
c	page 40-42 of MAST-811 Lecture Notes
c
c compiles with gfortran on OSX 10.11.6 (gcc version 5.4.0)
c	(MacBook Pro 2.5 GHz Intel Core i7)
c
c Andreas Muenchow, University of Delaware
c	Mar.-13, 2018
c
	integer n,Ntot,k
	real fc,dt,arg,pi,h(1000),hsum
	real hfreq,f,h2(1000),hfreq2
c
c Cut-Off Frequency
c
	fc = 0.3
	dt = 1.
	Ntot = 100
	pi = 4.*atan2(1.,1.)
	h0 = 2.*fc/dt
	write(16,*)'0 ',h0
	hsum = 0
	do 1 n=1,Ntot
c
c Truncated sinc function
c
	   arg = 2.*pi*float(n) * fc*dt
	   h(n) = 2.*fc*sin(arg)/arg
	   hsum = hsum+h(n)
c	   
c Lanczos corrections
c
	   arg2 = 2.*pi*float(n) / (2*Ntot+1)
	   h2(n) = sin(arg2)/arg2
c
c Write time-domain results to file for plotting
c
	   write(16,*)n,h(n),h2(n)*h(n)
1	continue
	write(6,*)hsum,h0+2*hsum
	do 2 k=1,Ntot
c
c Fourier Transform of truncated Sinc function
c
	   f = k/(Ntot*dt)
	   hfreq = h0
	   hfreq2 = h0
	   do 3 n=1,Ntot
		arg = 2.*pi*f*n*dt
		hfreq  = hfreq+2*h(n)*cos(arg)
c
c with Lanczos corrections
c
		hfreq2 = hfreq2+2*h(n)*cos(arg)*h2(n)
3	   continue
c
c Write frequency-domain results to file for plotting
c
		write(17,*)f,hfreq,hfreq2
2	continue
	stop
	end