c   MFPLOT
c
c   plots the current radially averaged log power spectrum with the
c   current model power spectrum.  Also plots the current bandpass
c   filters.
c
c   This is an auxilliary program in the matched filtering sequence.
c   It is intended to be used after mfinit and mfdesign.
c
c   input:  mfilt.9 radially averaged log power spectrum
c           mfilt.7 filter parameters
c
c   link with PLOT.LIB and GRAPHICS.LIB
c
c   Jeff Phillips July 1992
c
c      subroutine plot(selanar, iplotr)
      dimension rad(256), filt(256, 10), zp(10), bp(10), ilam(10),
     1 bound(10)
c
c  --  input results from mfinit
c
c      logical*1 selanar
      data pi / 3.14159265 /
      gmax = -100.
      gmin = 100.
      iplotr= ask_int('Enter plot device code (5=HP,8=CGA,9=EGA,10=VGA)
     &: ')
c      rewind(unit=9)
      open(9,file='mfilt.9',status='old',form='formatted',err=990)
      do 10 i = 1, 256
c      read(unit=9, fmt=*, err=990, end=11) rad(i), filt(i,1)
      read(unit=9, fmt=*, err=11, end=11) rad(i), filt(i,1)
      if (rad(i) .eq. 0) read(unit=9, fmt=*, end=11) rad(i), filt(i,1)
      if (filt(i,1) .gt. gmax) gmax = filt(i,1)
      if (filt(i,1) .lt. gmin) gmin = filt(i,1)
   10 continue
c
c  --  input filter parameters
c
   11 npts = i - 1
      open(7,file='mfilt.7',status='old',form='formatted',err=991)
      do 32 i = 1, 10
c      read(unit=7, fmt=*, err=991, end=33) zp(i), bp(i), ilam(i)
      read(unit=7, fmt=*, err=33, end=33) zp(i), bp(i), ilam(i)
   32 continue
   33 k = i - 1
      do 95 i = 1, npts
      filt(i,2) = 0.
      do 94 j = 1, k
      br = bp(j)
      if (ilam(j) .eq. 0) br = br * rad(i)
      if (ilam(j) .eq. 2) br = br / rad(i)
      fac = - (zp(j) * rad(i))
      if (fac .lt. (-88.028)) fac = -88.028
   94 filt(i,2) = filt(i,2) + (br * exp(fac))
   95 filt(i,2) = 2. * alog(abs(filt(i,2)))
c      call sendesc('[2J')
      call rplot(0., rad(npts), gmin, gmax, 2, iplotr)
      call curv(rad, filt(1,2), npts, 300)
      call curv(rad, filt(1,1), npts, 200)
      call endit
c      call sendesc(char(12))
c      if (selanar) call sendesc(char(50))
      do 96 j = 1, k
      write(unit=6, fmt=105) zp(j), bp(j), ilam(j)
   96 continue
  105 format(40x,3hzp=,f10.4,4h bp=,g10.4,6h ilam=,i3)
      write(unit=6, fmt=*) char(7)
      read(unit=5, fmt=*) 
      do 40 l = 1, k
      do 53 i = 1, npts
      filt(i,l) = 1.0
      if (ilam(l) - 1) 81, 82, 83
   81 bs = bp(l) * rad(i)
      goto 84
   82 bs = bp(l)
      goto 84
   83 bs = bp(l) / rad(i)
   84 continue
      do 53 il = 1, k
      if (il .eq. l) goto 53
      if (ilam(il) - 1) 101, 102, 103
  101 br = (bp(il) * rad(i)) / bs
      goto 104
  102 br = bp(il) / bs
      goto 104
  103 br = (bp(il) / rad(i)) / bs
c     if(arg*alog(abs(br)).gt.88.028) arg=88.028
  104 arg = (zp(l) - zp(il)) * rad(i)
      if (arg .gt. 88.028) then
      arg = 88.028
      br = amin1(br,1 / br)
      end if
      filt(i,l) = filt(i,l) + (br * exp(arg))
   53 continue
      do 40 i = 1, npts
      filt(i,l) = 1 / filt(i,l)
   40 continue
c      call sendesc('[2J')
      call rplot(0., rad(npts), 0., 1., 1, iplotr)
      do 50 i = 1, l
      call curv(rad, filt(1,i), npts, i*100)
   50 continue
      call endit
c      call sendesc(char(12))
c      if (selanar) call sendesc(char(50))
      close(unit=7)
      bound(1)=2.*pi/rad(npts) 		      
      do 54 j=1,l-1   
      bound(j+1)=0.
      do 51 i=1,npts+1
      if(filt(npts-i,j).gt.filt(npts-i,j+1)) go to 51
      bound(j+1)=2.*pi/rad(npts-i)
      go to 52
   51 continue
   52 write(unit=6, fmt=106) j, bound(j), bound(j+1)
   54 continue	     
c      bound(j)=2.*pi/rad(1)
c      print *, bound(j), bound(j+1) 
c      write(unit=6, fmt=106) (bound(i),i=1,l-1)
  106 format(1x,'bandpass wavelengths for layer ',i2,': '2(g10.4,3x))
      return 
  990 write(unit=*, fmt=*) 
     &' File not found:  use initialize and design before plot.'
      return 
  991 write(unit=*, fmt=*) ' File not found:  use design before plot.'
      end
      subroutine rplot(rmin, rmax, hmin, hmax, ityp, iplotr)
      dimension xp(4), yp(4), dxp(2), dyp(2)
      dimension x(n), y(n)
c      integer str(n)
      call pltset(iplotr, xp(4), yp(4), 1)
      xp(4) = min(xp(4),10.)
      yp(4) = min(yp(4),8.)
      xp(4) = xp(4) - .25
      yp(4) = yp(4) - .25
      dxp(1) = rmin
      dxp(2) = (rmax + .01) - mod(rmax,.01)
      dyp(1) = hmin
      dyp(2) = hmax
      xp(1) = xp(4) - 1.5
      xp(2) = 0.
      xp(3) = 1.25
      yp(1) = yp(4) - 1.25
      yp(2) = 0.
      yp(3) = 1.
      call scale(dxp, dyp, xp, yp, 4, ier)
      if (ier .ne. 0) goto 90
      goto (51, 52, 53), ityp
   51 call yaxis(dyp, dxp, yp, 0.1, 2, .12, '(f4.1)', 4)
      goto 54
   52 call yaxis(dyp, dxp, yp, 1.0, 2, .12, '(f4.0)', 4)
      goto 54
   53 call yaxis(dyp, dxp, yp, 0.5, 2, .12, '(f4.1)', 4)
   54 goto (55, 55, 56), ityp
   55 call xaxis(dxp, dyp, xp, 0.1, 2, .12, '(f4.1)', 4)
      goto 57
   56 call xaxis(dxp, dyp, xp, 5.0, 2, .12, '(f4.0)', 4)
   57 call neatl
      goto (58, 58, 59), ityp
   58 call vchar((dxp(1) + dxp(2)) * .5, dyp(1), 'radians/km', 10, 2, 
     &.12, 0., -.54, -.54)
      goto 60
   59 call vchar((dxp(1) + dxp(2)) * .5, dyp(1), 'azimuth', 7, 2, .12, 
     &0., -.54, -.54)
   60 goto (61, 62, 62), ityp
   61 call vchar(dxp, (dyp(1) + dyp(2)) * .5, 'amplitude', 9, 2, .12, 
     &1.5706, -.48, .54)
      goto 64
   62 call vchar(dxp, (dyp(1) + dyp(2)) * .5, 'log power', 9, 2, .12, 
     &1.5706, -.48, .54)
   64 return 
   90 stop 
      entry curv(x, y, n, m)
      call line(x, y, n, 0, m)
c      call line(x, y, 1, 0, m)
c      do i=2,n-1
c      call line(x(i),y(i),2,1,m)
c      end do
      return 
c      entry text(str,n,off)
c      call vchar((dxp(1) + dxp(2)) * .1, dyp(2) - off, str, n, 2, .10,
c     & 0., 0., 0.)
c      call home
c      return      
      entry endit()
      call endpt(ie)
      return 
      end
      function ask_int(request)
      character request*(*)
    2 write(unit=6, fmt=100) request
  100 format(1x,a,$)
      read(unit=5, fmt=*, err=1) ask_int
      return 
    1 write(unit=6, fmt=101) 
  101 format(/,20h **ERROR - try again)
      call bell
      goto 2
c-----------------------------------------------------------------------
c---
      end
      subroutine bell()
      character ding*1
      ding = char(7)
      write(unit=*, fmt=*) ding
      return 
c-----------------------------------------------------------------------
      end

