c....this one will make sure that ratio between bandwidth and frequency
c....is not below a certain value


      implicit double precision (a-h, o-z)
      double precision y(37),rootr(36),rooti(36)
      double complex zero,one
      real * 4 a(40)
      character * 16 anal,unstab
      double precision r(37),th(37)
      double complex ww, w(37)
      zero = (0.d0,0.d0)
      one = (1.d0,0.d0)
      write(6,110)
110   format(' Enter name of file containing analysis data.')
      read(5,11) anal
11    format(a16)
      write(6,112)
112   format(' Enter range of frames to correct.') 
      read(5,610) ifirst,ilast
610   format(2i6)
      write(6,113)
113   format(' Enter number of poles (max is 36).')
      read(5,114) npoles
114   format(i2)
      write(6,600)
600   format(' Enter sampling rate.')
      read(5,601) sr
601   format(f8.0)
      write(6,602)
602   format( ' Enter maximum ration of bw to freq')
      read(5,603) bandmx
603   format(f8.0)
      ndata = 4
      nlocs = (npoles + ndata) * 4
      pi = 3.141592654
      fynq = .5 * sr / pi
      bwfact = sr/pi
      print*, npoles,nlocs,sr,fynq,bwfact


	open(unit=18, file=anal,status='old',
     &	access='direct',form='unformatted',recl=nlocs)
      k4=npoles+1
      k4m=k4-1
      nall = k4m + ndata
      do 199 kkk = ifirst,ilast
      read(18,rec=kkk,end=999) (a(ll),ll=1,nall)
      write(6,198) kkk
198   format( ' correcting frame  ',i5)
      do 1601 ii=1,k4m
c     print*, a(ii+4)
1601  y(ii)= -a(ii+4)
      y(k4)=1.
c       find out if we want to suppress printing
      kprint=0
      eps=10.d0**(-8)
303   call factor(y,k4,rootr,rooti,kinsid,kprint,eps)
c     write(6,916)kinsid
916   format('there are',i5,'roots inside the unit circle')
      do 100 j=1,k4m
      r(j) = dsqrt(rootr(j) **2 + rooti(j)**2)
      th(j) = datan2(rooti(j),rootr(j))
      freq = th(j)*fynq
      if(freq.eq.0.) goto 100
      bw = (1.d0 - r(j))*bwfact
      ratio = bw/abs(freq)
c      print*,bw,freq,ratio,r(j),j
      if (ratio.gt.bandmx) goto 100
      xbw = bandmx * abs(freq)
      r(j) = 1.0 - xbw/bwfact
c      print*,xbw,r(j)
100   continue 
      do 10 k=1,k4m
10    w(k) = zero
      w(k4)=one
      do 20 k=1,k4m
c	ww=dcmplx(rootr(k),rooti(k))
	ww=dcmplx(r(k)*dcos(th(k)),r(k)*dsin(th(k)))
	l1=k4-k
	do 12 j=l1,k4m
12    	   w(j)=w(j+1)-ww*w(j)
20      w(k4)=-ww*w(k4)
	do 30 j=2,k4 
	zz=real(w(j))
	a(k4+ndata+1-j) = -zz
30	continue
        write(18,rec=kkk) (a(ll),ll=1,nall)
199     continue
999    stop 
      end
