c
c Illustrations for Lanczos Filter as discussed
c	page 40-42 of MAST-811 Lecture Notes
c
c Andreas Muenchow, University of Delaware
c	Sept.-26, 2019
c	Mar.-13, 2018
c
	integer n,Ntot,k
	real fc,dt,arg,pi,h(2000),hsum,hsum2
	real hfreq,f,h2(2000),hfreq2
c
c Cut-Off Frequency
c
	fc = 0.3
	dt = 1.
	Ntot = 48
	write(6,*)'Enter Half-Window-Width (hrs)'
	read(5,*)Ntot
	pi = 4.*atan2(1.,1.)
	h0 = 2.*fc/dt
	write(16,*)'0 ',h0
	hsum = 0.
	hsum2 = 0.
	do 1 n=1,Ntot
c
c Truncated sinc function
c
	   arg  = 2.*pi*fc*float(n)*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
	   hsum2 = hsum2+h2(n)*h(n)
c
c Write time-domain results to file for plotting
c
	   write(16,*)n,h(n),h2(n)*h(n)
1	continue
	write(6,*)'Sum of Filter Weights without Lanczos corrections: ',h0+2*hsum
	write(6,*)'Sum of Filter Weights with    Lanczos corrections: ',h0+2*hsum2
c
c Calculate Fourier Transform of Filter Weights for "contineous" frequencies
c
	do 2 k=1,1001
c	do 2 k=1,Ntot
c
c Fourier Transform of truncated Sinc function
c
	   f = float(k-1)/(float(Ntot)*dt)
	   f = f*0.1 
	   f = (k-1)*0.5/1000.
	   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