c
c Development of Lanczos Raised Cosine Filter
c	Sept.-30, 2010 amue
c
c There is a problem with a small offset at f<fc that depends
c	on the number of filter weights when it should oscillate
c	around 1. (<0.8% for s=100). Needs further checking and
c	verification.
c
c
	parameter (nn=2000)
	integer n,s
	double precision h(nn),ffilter(nn),ffilter2(nn),arg
	double precision hsum,htot
	double precision fc,dt,t
	pi = 4.*atan2(1.,1.)
	write(6,*)'Number of filter weights (<=500)'
	read(5,*)nweights
c
c Cut-off Period (10 hours)
c
99	continue
	t = 10.
	write(6,*)'Enter cut-off period (<50 hrs)'
	read(5,*)t
	if (t.gt.50) goto 99
c
c Cut-Off Frequency (0.1 cph)
c
	fc = 1./T
c
c Time step (1 hr)
c
	dt = 1.
c
c Filter Length
c
	s = nweights
c
c generate filter-weights from sinc function
c
	h(1) = 2.*fc/dt
	hsum = 0.
	do 1 n=2,s
	   arg = 2.*pi*fc*float(n-1)*dt
	   h(n) = 2.*fc*sin(arg)/arg
	   hsum = hsum+h(n)
1	continue
	htot = h(1)+2.*hsum
c
c Normalize the filter weights so they add up to 1
c
	h(1) = h(1)/htot
	hsum = 0.
	do 4 n=2,s
	   h(n) = h(n)/htot
	   hsum = hsum+h(n)
4	continue
	write(6,*)h(1),htot,h(1)+2.*hsum
	write(16,101)(n-1,h(n),n=1,s)
c
c Fourier Transform of filter-weights
c
c   Calculate for 2000 frequencies (just for plotting, no physical significance)
c
	fstep = 0.0001
	do 2 j=1,2000
c
c Raw filter to approximate ideal filter (composed of cosines)
c
	   ffilter(j) = h(1)
c
c Smoothed filter to better approximate filter (smoothed in frequency)
c
	   ffilter2(j) = h(1)
c
c Frequency (arbitrary, can calculate this for all frequencies
c
	   f = float(j-1)*fstep
	   do 3 n=2,s
	      arg =  2.*pi*float(n-1)*f*dt
	      arg2 = 2.*pi*float(n-1)/(2.*float(s-1)+1.)
	      ffilter(j)  = ffilter(j) +2.*h(n)*cos(arg)
	      ffilter2(j) = ffilter2(j)+2.*h(n)*cos(arg)*sin(arg2)/arg2
3	   continue
	   write(17,102)f,ffilter(j),ffilter2(j)
2	continue
101	format(i5,f20.8)
102	format(3f20.8)
	stop
	end