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       Oct.-7, 2021
c	Sept.-26, 2019
c	Mar.-13, 2018
c
	integer n,Ntot,k
	real fc,dt,arg,pi,h(2000),hsum,hsum2,h0
	real hfreq,f,h2(2000),hfreq2
	real h3(2000),hfreq3,h03
	real h4(2000),hfreq4,h04
c        
        idev = 1
        write(6,*)'Enter: 1 no correction, 0 with correction'
        read(5,*)idev
c
c Cut-Off Frequency
c
	fc = 0.1
        write(6,*)'Enter cut-off frequency in [0.1,0.4]'
        read(5,*)fc
	dt = 1.
	Ntot = 48
	write(6,*)'Enter Half-Window-Width (hrs)'
	read(5,*)Ntot
	pi = 4.*atan2(1.,1.)
	h0 = 2.*fc/dt
        h03 = 1./float(Ntot)
        h04 = 1./float(2*Ntot+1)
	write(6,*)'0 ',h0,h0,h03
	write(16,*)'0 ',h0,h0,h03
	hsum = 0.
	hsum2 = 0.
	hsum3 = 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) / float(2*Ntot+1)
	   h2(n) = sin(arg2)/arg2
	   hsum2 = hsum2+h2(n)*h(n)
c
        if (idev.eq.1) then
           h2(n) = 1.
        end if

c
        h3(n) = ((1.+cos(2.*pi*n/(2*Ntot+1.))))/(2.*Ntot)
	   hsum3 = hsum3+h3(n)
c
           h4(n) = 1./float(2*Ntot+1)
	   hsum4 = hsum4+h4(n)
c
c Write time-domain results to file for plotting
c
	   write(16,*)n,h(n),h2(n)*h(n),h3(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
	write(6,*)'Sum of Filter Weights         Raised Cosine:       ',h03+2*hsum3
	write(6,*)'Sum of Filter Weights         Box Car: ',h04+2*hsum4
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
	   hfreq3 = h03
	   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)
                hfreq3 = hfreq3+2.*h3(n)*cos(arg)
3	   continue
c
c Write frequency-domain results to file for plotting
c
c               if (k .eq. 1) then
c                  hfreq3 = 1.
c               else
c                  hfreq3 = sin(arg)/(arg)
c               end if
		write(17,*)f,hfreq,hfreq2,hfreq3
		write(18,*)Ntot/10*f/fc,hfreq,hfreq1,hfreq3
2	continue
	stop
	end