c   MFDESIGN
c
c   interactive design of matched bandpass filters
c
c   input:  for009 log averaged radial power spectrum
c   output: for007 filter parameters
c
c   the following subroutines must be modified if array dimensions
c   are changed:  nlctl, parshl, marq, corvec
c
c      subroutine design(writeover, selanar, iplotr)
      common /work1/ npts, k, zp(10), bp(10), ilam(10), rad(256)
      dimension g(256), cs(256), ct(256), cu(256), bound(10)
      dimension filt(256, 10)
      character ans*1, ask_yn*1
c      logical writeover, selanar
      data dval / 1.e+38 /, pi/3.14159265/
      iplotr= ask_int('Enter plot device code (5=HP,8=CGA,9=EGA,10=VGA)
     &: ')
    3 gmin = 100.
c
c  --  input results from program mspec
c	   g(i) is the log averaged radial power spectrum (from init)
c	   cs(i) is the log residual radial power spectrum (effects of 
c	         previous layers removed)
c	   ct(i) is the log reduced residual radial power spectrum (corrected
c	         for current layer type)
c	   cu(i) is the modeled log radial power spectrum for the current
c	         layer
c
c	it=1
      gmax = -100.
c      rewind(unit=9)
      open(9,file='mfilt.9',form='formatted',status='old',err=990)
      do 10 i = 1, 256
c      read(unit=9, fmt=*, err=990, end=11) rad(i), g(i)
      read(unit=9, fmt=*, err=11, end=11) rad(i), g(i)
      if (rad(i) .eq. 0.) read(unit=9, fmt=*, end=11) rad(i), g(i)
      if (g(i) .gt. gmax) gmax = g(i)
      if (g(i) .lt. gmin) gmin = g(i)
      cs(i) = g(i)
      ct(i) = g(i)
   10 continue
   11 npts = i - 1
      gmax = float(int(gmax + 1.))
      gmin = float(int(gmin - 1.))
      hmax = gmax
      hmin = gmin
      fmax = gmax
c	if(it.ne.1) go to 93
      fmin = gmin
      k = 1
      jcont = -99
c
c  --  enter layer code
c
      icont = 0
    2 if (icont .ne. (-1)) write(unit=6, fmt=100) 
  100 format('  0  magnetic dipole layer',/,
     &'  1  magnetic half-space or density layer',/,
     &'  2  density half-space',/,' -1  next layer',/,' 99  exit',/,/)
      write(unit=6, fmt=101) 
  101 format(' enter layer code: '$)
      read(unit=5, fmt=*, err=2) icont
      if (icont .eq. jcont) goto 15
      jcont = icont
      fmax = hmax
      fmin = hmin
      if ((icont .lt. 0) .or. (icont .gt. 2)) goto 91
   12 fmax = -100.
c
c  --  perform layer correction
c
      fmin = 100.
      sign = 1 - icont
      do 13 i = 1, npts
      ct(i) = cs(i) - ((sign * 2.) * alog(rad(i)))
      if (ct(i) .gt. fmax) fmax = ct(i)
      if (ct(i) .lt. fmin) fmin = ct(i)
   13 continue
      fmax = float(int(fmax + 1.))
      fmin = float(int(fmin - 1.))
   14 continue
      ilam(k) = icont
c   15 call sendesc('[2J')
   15 continue
      call rplot(0., rad(npts), fmin, fmax, 2, iplotr)
      call curv(rad, ct, npts, 300)
c      if (.not. writeover) then
c      write(unit=6, fmt=*) char(7)
c      read(unit=5, fmt=*)
c
c  --  input parameters for current layer
c
c      end if
c      call sendesc('[H')
      write(unit=6, fmt=*) 
     &' input: left-intercept, bottom-intercept (or 0)'
c      call text('input: left-intercept, bottom-intercept (or 0)',46,2.0)
      read(unit=5, fmt=*, err=15) br, h1
      h0 = fmin
      if (h1 .gt. 0) goto 16
      write(unit=6, fmt=*) ' input: right-intercept'
c      call text('input: right-intercept',22,3.0)
      read(unit=5, fmt=*, err=15) h0
      h1 = rad(npts)
   16 h0 = br - h0
      zp(k) = abs((h0 * .5) / h1)
      bp(k) = exp(br * .5)
      do 30 i = 1, npts
      cu(i) = 0.
      fac = - (zp(k) * rad(i))
      if (fac .lt. (-88.028)) fac = -88.028
      br = bp(k)
      fac = abs(br * exp(fac))
      if (fac .lt. 5.888e-39) fac = 5.888e-39
      cu(i) = 2. * alog(fac)
c
c  --  plot the modeled spectrum
c
   30 continue
      call curv(rad, cu, npts, 200)
      call endit
c      call sendesc(char(12))
c      if (selanar) call sendesc(char(50))
      write(unit=6, fmt=102) zp(k), bp(k), ilam(k)
  102 format(25x,'depth=',g10.4,' amplitude=',g10.4,' layer code=',i3)
      goto 2
   91 continue
c
c -- compute residual power spectrum (remove the effects of the previous
c layer)
c
      if ((k .eq. 10) .or. (icont .gt. 2)) goto 93
      hmax = -100
      hmin = 100
      do 92 i = 1, npts
      fac = - (zp(k) * rad(i))
      if (fac .lt. (-88.028)) fac = -88.028
      br = bp(k)
      if (ilam(k) .eq. 0) br = br * rad(i)
      if (ilam(k) .eq. 2) br = br / rad(i)
      cs(i) = 2. * alog(abs(exp(cs(i) * .5) - (br * exp(fac))))
      ct(i) = cs(i)
      if (cs(i) .gt. hmax) hmax = cs(i)
      if (cs(i) .lt. hmin) hmin = cs(i)
   92 continue
      hmax = float(int(hmax + 1.))
      hmin = float(int(hmin - 1.))
      fmax = hmax
      fmin = hmin
      k = k + 1
      goto 2
c	it=2
c
c -- plot the observed and modeled specta
c
   93 continue
      call compmod(npts, k, bp, zp, ilam, rad, cs)
c      call sendesc('[2J')
      call rplot(0., rad(npts), gmin, gmax, 2, iplotr)
      call curv(rad, cs, npts, 300)
      call curv(rad, g, 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=102) zp(j), bp(j), ilam(j)
c	call askyn('Do you want to iterate on this solution?: ',ans)
   96 continue
      ans = ask_yn('Do you want to iterate on this solution?: ')
      if (ans .eq. 'n') goto 97
      call nlctl(g, cs)
c97	call askyn('Do you want to start over?: ',ans)
      goto 93
   97 ans = ask_yn('Do you want to start over?: ')
c		it=1
      if (ans .eq. 'y') then
      rewind(unit=9)
      goto 3
      end if
      open(7,file='mfilt.7',form='formatted',status='unknown')
      do 98 j = 1, k
   98 write(unit=7, fmt=*) zp(j), bp(j), ilam(j)
      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) 201, 202, 203
  201 br = (bp(il) * rad(i)) / bs
      goto 204
  202 br = bp(il) / bs
      goto 204
  203 br = (bp(il) / rad(i)) / bs
  204 arg = (zp(l) - zp(il)) * rad(i)
      if (arg .gt. 88.028) then
      arg = 88.028
      br = amin1(br,1 / br)
      end if
      arg = br * exp(arg)
      if (arg .gt. (dval - filt(i,l))) goto 53
      filt(i,l) = filt(i,l) + arg
   53 continue
      do 40 i = 1, npts
c      print *,i, l, filt(i,l)
      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, k
      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))
   99 stop
  990 write(unit=*, fmt=*) 
     &'File not found: run mfinit before mfdesign.'
      stop
      end
      subroutine compmod(npts, k, bp, zp, ilam, rad, cs)
      dimension bp(k), zp(k), ilam(k), rad(npts), cs(npts)
      do 95 i = 1, npts
      cs(i) = 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 cs(i) = cs(i) + (br * exp(fac))
   95 cs(i) = 2. * alog(abs(cs(i)))
      return 
c***********************************************************************
c*******
      end
c  improve the model using non-linear least squares
      subroutine nlctl(g, cs)
      common /work1/ npts, k, zp(10), bp(10), ilam(10), rad(256)
      common /swork/ a(2560), err(256), s(20), u(2560), v(256)
      common /free1/ rms1, iter, xk1
      dimension g(npts), cs(npts)
      character ans*1
      logical step
      data zero / 0.0 /
      deltz = zp(1) / 10.
c
      deltb = bp(1) / 10.
      write(unit=*, fmt=175) 
  175 format(35h step-by-step or forever (s or f)? ,$)
      read(unit=5, fmt=103) ans
c  get error vector
      if ((ans .eq. 's') .or. (ans .eq. 'S')) step = .true.
      iter = 1
  100 call rmserr(npts, rms1, g, cs, err)
      if (step) write(unit=*, fmt=176) rms1
c  setup jacobian
  176 format(21h  starting rms error ,1pg15.5)
      iad = 1
      do 20 i = 1, (npts * 2) * k
   20 a(i) = 0.0
      do 30 i = 1, k
      call parshl(deltb, zero, i, a(iad))
   30 iad = iad + npts
      do 40 i = 1, k
      call parshl(zero, deltz, i, a(iad))
c  scale jacobian
   40 iad = iad + npts
      call cscale(npts, 2 * k, a, s)
      if (iter .gt. 1) goto 200
      if (step) call crlmtx(6, npts, 2 * k, a)
  200 call marq(step, g, *300)
      iter = iter + 1
      if (step) then
      write(unit=*, fmt=210) 
  210 format(20h another iteration? $)
      read(unit=5, fmt=103) ans
  103 format(a1)
      if ((ans .eq. 'n') .or. (ans .eq. 'N')) goto 300
      end if
      call compmod(npts, k, bp, zp, ilam, rad, cs)
      goto 100
  300 return 
c
      end
c  partial derivatives via central differences
      subroutine parshl(db, dz, iparm, pf)
      common /work1/ npts, k, zp(10), bp(10), ilam(10), rad(256)
      dimension tmp(256), pf(256)
      dimension bp1(10), zp1(10)
      if ((iparm .lt. 1) .or. (iparm .gt. k)) goto 400
      do 10 i = 1, npts
      pf(i) = 0.0
   10 tmp(i) = 0.0
      if (db .eq. 0.0) goto 200
      do 20 i = 1, k
   20 bp1(i) = bp(i)
      bp1(iparm) = bp1(iparm) + db
      call compmod(npts, k, bp1, zp, ilam, rad, pf)
      bp1(iparm) = bp1(iparm) - (2. * db)
      call compmod(npts, k, bp1, zp, ilam, rad, tmp)
      delta = 1.0 / (2.0 * db)
      goto 300
  200 if (dz .eq. 0.0) goto 400
      do 30 i = 1, k
   30 zp1(i) = zp(i)
      zp1(iparm) = zp1(iparm) + dz
      call compmod(npts, k, bp, zp1, ilam, rad, pf)
      zp1(iparm) = zp1(iparm) - (2. * dz)
      call compmod(npts, k, bp, zp1, ilam, rad, tmp)
      delta = 1.0 / (2.0 * dz)
  300 do 40 i = 1, npts
   40 pf(i) = delta * (pf(i) - tmp(i))
      return 
  400 write(unit=*, fmt=*) 'Error in parshl parameters'
      return 
      end
c  print the triangular correlation matrix
      subroutine crlmtx(jdev, m, n, a)
      dimension icl(20), a(m, n)
      write(unit=jdev, fmt=5) 
    5 format(1x,3h a=)
      do 20 j = 2, n
      do 10 i = 1, j - 1
   10 icl(i) = int((100.0 * correl(m,a(1,j),a(1,i))) + 0.5)
      write(unit=jdev, fmt=15) j, (icl(i),i = 1, j - 1)
   15 format(1x,i3,20i3)
c	write(jdev,30)
c30	format(/,1x,' b=  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15')
   20 continue
      return 
c***********************************************************************
c*******
      end
c  column scaling for matrix a
      subroutine cscale(m, n, a, s)
      dimension s(n), a(1)
      tol = 1.e-10
      do 10 j = 1, n
      k = ((j - 1) * m) + 1
      pmin = a(k)
      pmax = a(k)
      do 5 i = 1, m
      if (a(k) .gt. pmax) pmax = a(k)
      if (a(k) .lt. pmin) pmin = a(k)
    5 k = k + 1
      s(j) = amax1(abs(pmax),abs(pmin))
      if (s(j) .lt. tol) s(j) = 1.0
   10 s(j) = 1.0 / s(j)
      k = 1
      do 20 j = 1, n
      do 20 i = 1, m
      a(k) = a(k) * s(j)
   20 k = k + 1
      return 
      end
      subroutine marq(step, g, *)
      common /swork/ a(2560), err(256), s(20), u(2560), v(256)
      common /work1/ npts, k, zp(10), bp(10), ilam(10), rad(256)
      common /free1/ rms1, iter, xk1
c
      dimension g(256), pf(256), err1(256)
      dimension q(20), work(80), ute(20), qute(20), dp(20), dg(20)
      equivalence (ute, work), (qute, work(21)), (dp, work(41)), (dg, 
     &work(61))
      dimension bp1(10), zp1(10)
c
c     dimension nf1(10)
      dimension xk(3), rms(3), ang(3), dp2(3)
c	logical grvswt,magswt
      character cp*1
      logical step
c  eigenvalue attenuation rate set by attn
c
c	do 10 i=1,5
c	nf1(i)  = nf(i)  +(igs-1)
c10	nf1(i+5)= nf(i+5)+(ihs-1)
      data attn / 4. /
      data mxridg / 10 /
      np = 2 * k
      icount = 0
      istop = 0
      limit = 0
      init = 1
c
      ii = 3
      call svd(a, npts, np, q, u, v, work)
      call atb(u, npts, np, err, npts, 1, ute)
      if (step) write(unit=*, fmt=*) ' eigenvalues'
      if (step) write(unit=*, fmt=*) (q(i),i = 1, np)
      call ata(u, npts, np, a)
      call chki(a, np, icode)
      call ata(v, np, np, a)
      call chki(a, np, jcode)
      if (icode .lt. 6) write(unit=*, fmt=*) 'svd: UtU', icode
c
      if (jcode .lt. 6) write(unit=*, fmt=*) 'svd: VtV', jcode
      call ammi(np, q, qmin, qmax, mn, mx, 1)
      if ((init .eq. 1) .and. (ii .eq. 3)) xk(3) = qmax
c
c  steepest descent vector
      if (iter .gt. 1) xk(3) = xk1 * 1.2
      damp = 10. * (qmax ** attn)
      do 50 i = 1, np
   50 qute(i) = ute(i) / (q(i) + damp)
      call axb(v, np, np, qute, np, 1, dg)
      dgl = 0.
      do 80 i = 1, np
      dg(i) = dg(i) * s(i)
   80 dgl = dgl + (dg(i) * dg(i))
c
c  ridge regression loop
      dgl = sqrt(dgl)
      if (step) write(unit=*, fmt=90) 
   90 format(31h   cutoff	  rms	angle    ave dp)
  100 icount = icount + 1
      if (((limit .eq. 1) .and. (init .eq. 0)) .or. (icount .gt. mxridg)
     &) then
      istop = 1
      call ammi(3, rms, emin, emax, ii, mx, 1)
      call ammi(3, ang, amin, amax, iia, mx, 1)
      if (iia .gt. ii) ii = iia
      if (step) write(unit=*, fmt=101) xk(ii)
  101 format(27h choosing eigenvalue cutoff,g15.4)
c
c  get new correction vector
c	if(manual.eq.1) then
c	  type *,'damping ='
c	  read(5,*) xk(ii)
c	endif
      end if
      damp = xk(ii) ** attn
      do 104 i = 1, np
      q2 = q(i) ** (attn - 1.)
  104 qute(i) = ute(i) * (q2 / ((q2 * q(i)) + damp))
      call axb(v, np, np, qute, np, 1, dp)
      do 105 i = 1, np
c
c  assemble and check parameter vector
  105 dp(i) = dp(i) * s(i)
c
c  forward calculation
      call corvec(limit, step, dp, bp1, zp1)
  200 call compmod(npts, k, bp1, zp1, ilam, rad, pf)
c	if(grvswt) call mardvr(0,xb1,zb1,rho1,igs,ngrv,
c	1	    f(nf1(4)),f(nf1(5)),rmsg)
c	if(magswt) call mardvr(1,xb1,zb1,sus1,ihs,nmag,
c	1	    f(nf1(9)),f(nf1(10)),rmsm)
      call rmserr(npts, rmsg, g, pf, err1)
c	if(magswt) rms(ii)=rmsm
c	if(grvswt .and. magswt) rms(ii)=(weight*rmsg+rmsm)/(weight+1.0)
      rms(ii) = rmsg
      dp2(ii) = 0.0
      an = 0.0
      dpl = 0.0
      do 210 i = 1, np
      dp2(ii) = dp2(ii) + (dp(i) * dp(i))
      dpl = dpl + (dp(i) * dp(i))
  210 an = an + (dp(i) * dg(i))
      dpl = sqrt(dpl)
      if (dpl .lt. 1.e-10) goto 300
      ad = dpl * dgl
      ang(ii) = acosd2(an,ad)
      dp2(ii) = sqrt(dp2(ii) / float(np))
      if (step) write(unit=*, fmt=220) xk(ii), rms(ii), ang(ii), dp2(ii)
  220 format(1x,2g12.4,f7.2,f10.3)
c
c  iteration control
      if (istop .eq. 1) goto 300
      if ((init .eq. 1) .and. (ii .gt. 1)) then
      ii = ii - 1
      xk(ii) = xk(ii + 1) / 2.
      else
      init = 0
      call bictl(xk, rms, ang, dp2, ii, ireset)
      if (ireset .eq. (-1)) limit = 1
      end if
c
c  print new parameter values
      goto 100
  300 if (step) write(unit=*, fmt=310) 
  310 format(25h parm   change	 new value)
      do 390 i = 1, np
      if (i .gt. k) goto 320
      cp = 'bp'
      j = i
      p = bp1(i)
      goto 380
  320 cp = 'zp'
      j = i - k
      p = zp1(i - k)
  380 if (step) write(unit=*, fmt=381) cp, j, dp(i), p
  381 format(1x,a2,i2,3x,1p2(g12.4,1x))
c
c  save new parameter values
  390 continue
      rmsimp = 0.0
      if (rms1 .gt. 0.0) rmsimp = (100. * (rms1 - rms(ii))) / rms1
      write(unit=*, fmt=410) iter, rms(ii), rmsimp
  410 format(11h iteration ,i2,17h,  rms error is  ,1pg15.4,/,
     &25h percent improvement is  ,1pg15.4)
      if (step) then
      write(unit=*, fmt=420) 
  420 format(28h save new parameter values ?$)
      read(unit=5, fmt=103) ans
  103 format(a1)
      if ((ans .eq. 'n') .or. (ans .eq. 'N')) return 
      else
      if (rmsimp .le. 0.0) return 1
      end if
      do 430 i = 1, k
      bp(i) = bp1(i)
  430 zp(i) = zp1(i)
      xk1 = xk(ii)
      rms1 = rms(ii)
      return 
c
      end
      subroutine ammi(n, a, amn, amx, mn, mx, isave)
      dimension a(n)
      if (isave .gt. 1) goto 1
      mn = 1
      mx = 1
      amn = a(1)
      amx = a(1)
    1 do 3 i = 1, n
      if (a(i) .ge. amn) goto 2
      amn = a(i)
      mn = i
    2 if (a(i) .le. amx) goto 3
      amx = a(i)
      mx = i
    3 continue
      return 
c-----------------------------------------------------------------------
c-
      end
      subroutine ata(a, m, n, b)
      dimension a(m, n), b(n, n)
      do 10 j = 1, n
      do 10 i = 1, n
   10 b(i,j) = xty(a(1,i),1,a(1,j),1,m)
      return 
      end
      subroutine atb(a, ma, na, b, mb, nb, c)
      dimension a(ma, na), b(mb, nb), c(na, nb)
      if (ma .ne. mb) goto 20
      k = ma
      do 10 j = 1, nb
      do 10 i = 1, na
   10 c(i,j) = xty(a(1,i),1,b(1,j),1,k)
      return 
   20 write(unit=*, fmt=21) 
   21 format(16h atb: col.ne.row)
      return 
      end
      subroutine axb(a, ma, na, b, mb, nb, c)
      dimension a(ma, na), b(mb, nb), c(ma, nb)
      if (na .ne. mb) goto 20
      k = na
      do 10 j = 1, nb
      do 10 i = 1, ma
   10 c(i,j) = xty(a(i,1),ma,b(1,j),1,k)
      return 
   20 write(unit=*, fmt=21) 
   21 format(22h axb: col a .ne. row b)
      return 
c***********************************************************************
c*******
      end
c  locate the next eigenvalue cutoff by checking the rms error and the a
cngle
c  formed between test vectors and the steepest descent vector and choos
ce
c  the more conservative minimum.
c  ireset=-1 stop,		 ireset=0 interval division,
c  ireset=1 interval addition, ireset=2 balance intervals
      subroutine bictl(x, err, ang, dp2, ii, ireset)
      dimension x(3), err(3), ang(3), dp2(3)
      ireset = -1
c  check for convergence
      if ((x(3) .le. x(2)) .or. (x(2) .le. x(1))) return 
      if (err(2) .eq. 0.0) return 
      derr = (abs(err(3) - err(2)) + abs(err(2) - err(1))) / err(2)
c  balance intervals
      if (derr .lt. 0.001) return 
      ratio = (x(3) - x(2)) / (x(2) - x(1))
      if ((ratio .gt. 5.0) .or. (ratio .lt. 0.2)) then
      ireset = 2
      ii = 2
      x(2) = 0.5 * (x(3) + x(1))
      goto 100
c
      end if
      icase = 0
      se2 = (err(2) - err(1)) / (x(2) - x(1))
      se3 = (err(3) - err(2)) / (x(3) - x(2))
      if ((se2 .gt. 0.) .and. (se3 .gt. 0.0)) icase = 1
      if ((se2 .lt. 0.) .and. (se3 .lt. 0.0)) icase = 4
      if (icase .eq. 0) then
      icase = 2
      if (abs(se2) .gt. se3) icase = 3
c
      end if
      jcase = 0
      sa2 = (ang(2) - ang(1)) / (x(2) - x(1))
      sa3 = (ang(3) - ang(2)) / (x(3) - x(2))
      if ((sa2 .gt. 0.) .and. (sa3 .gt. 0.0)) jcase = 1
      if ((sa2 .lt. 0.) .and. (sa3 .lt. 0.0)) jcase = 4
      if (jcase .eq. 0) then
      jcase = 2
      if (abs(sa2) .gt. sa3) jcase = 3
c
c	write(33,*) 
c	write(33,*) ' cutoff',x
c	write(33,*) ' err   ',err
c	write(33,*) ' ang   ',ang
c	write(33,*) ' i/j c ',icase,jcase
c
c  override an interval division by a good angle improvement
      end if
      dangl = abs(ang(3) - ang(2)) + abs(ang(2) - ang(1))
c
      if ((icase .le. 2) .and. ((jcase .gt. 2) .and. (dangl .gt. 5.0))) 
     &icase = (icase + jcase) / 2
      ireset = 1
      if (icase .eq. 1) then
      x(3) = x(2)
      x(2) = x(1)
      x(1) = x(2) / 2.
      err(3) = err(2)
      err(2) = err(1)
      ang(3) = ang(2)
      ang(2) = ang(1)
      dp2(3) = dp2(2)
      dp2(2) = dp2(1)
      else if ((icase .eq. 2) .or. (icase .eq. 3)) then
      if (icase .eq. 2) x(1) = (x(1) + x(2)) / 2.
      if (icase .eq. 3) x(3) = (x(2) + x(3)) / 2.
      ireset = 0
      else if (icase .eq. 4) then
      x(1) = x(2)
      x(2) = x(3)
      x(3) = (1.5 * (x(2) - x(1))) + x(2)
      err(1) = err(2)
      err(2) = err(3)
      ang(1) = ang(2)
      ang(2) = ang(3)
      dp2(1) = dp2(2)
      dp2(2) = dp2(3)
      end if
      ii = 1
      if (icase .gt. 2) ii = 3
  100 err(ii) = 1.0e30
c	write(33,*) ' bictl:',ii,x(ii)
      ang(ii) = 90.0
      return 
      end
c  check the identity matrix
      subroutine chki(a, n, iexp)
      dimension a(n, n)
      iexp = 15
      tmx = 0.0
      do 10 j = 1, n
      do 10 i = 1, n
      t = abs(a(i,j))
      if (i .eq. j) t = abs(1.0 - a(i,i))
      if (t .gt. tmx) tmx = t
   10 continue
      if (tmx .le. 1.e-15) return
      dum=alog10(tmx) 
      iexp = int(0.5 - dum)
      return 
c
      end
      subroutine corvec(limit, step, dp, bp1, zp1)
c
      common /work1/ npts, k, zp(10), bp(10), ilam(10), rad(256)
      dimension dp(20), bp1(10), zp1(10)
c	data rlimit/.05/	
c	logical grvswt,magswt,limitp
c
      logical step
c	x1=0.0
c	x2=0.0
c	if(grvswt) x1 = abs(f(nf(2)-1)-f(nf(1)) )/50.0
c	if(magswt) x2 = abs(f(nf(7)-1)-f(nf(6)) )/50.0	
c	xlimit=amax1(x1,x2)
c	z1=0.0
c	z2=0.0
c	if(grvswt) z1 = abs(f(nf(3)-1)-f(nf(2)) )/10.0
c	if(magswt) z2 = abs(f(nf(8)-1)-f(nf(7)) )/10.0	
c	zlimit=amax1(z1,z2)
      limit = 0
      blimit = abs(bp(k)) / 10.
c
c  assemble parameter vector
      zlimit = abs(zp(k) - zp(1)) / 10.
      do 10 i = 1, k
      bp1(i) = bp(i)
c
c  add in correction vector and check excursion
   10 zp1(i) = zp(i)
c  bp parameter
      ip = 0
      do i = 1, k
      limitb = 0
      ip = ip + 1
   50 absdp = abs(dp(ip))
      bp1(i) = bp(i) + dp(ip)
      if ((absdp .gt. blimit) .or. (bp1(i) .le. 0.0)) then
      dp(ip) = dp(ip) / 2.0
      if (absdp .lt. 0.01) dp(ip) = 0.0
      limitb = ip
      dpb = dp(ip)
      if (dpb .ne. 0.0) goto 50
      end if
      if (limitb .gt. 0) then
      if (step) write(unit=*, fmt=*) 
     &' limited excursion of bp parameter', limitb, ' to', dpb
      limit = 1
      end if
c  zp parameter
      end do
      do i = 1, k
      limitz = 0
      ip = ip + 1
  100 absdp = abs(dp(ip))
      zp1(i) = zp(i) + dp(ip)
      if ((absdp .gt. zlimit) .or. (zp1(i) .le. 0.0)) then
      dp(ip) = dp(ip) / 2.0
      if (absdp .lt. 0.01) dp(ip) = 0.0
      limitz = i
      dpz = dp(ip)
      if (dp(ip) .ne. 0.0) goto 100
      end if
      if (limitz .gt. 0) then
      if (step) write(unit=*, fmt=*) 
     &' limited excursion of zp parameter', limitz, ' to', dpz
      limit = 1
      end if
c
      end do
      return 
c***********************************************************************
c*******
      end
      subroutine rmserr(npts, errt, g, cs, err)
      dimension g(npts), cs(npts), err(npts)
      errt = 0.0
      do 10 i = 1, npts
      err(i) = g(i) - cs(i)
   10 errt = errt + (err(i) ** 2)
      errt = sqrt(errt / float(npts))
      return 
      end
      subroutine svd(a, m, n, s, u, v, wrk)
      dimension a(m, n), s(n), u(m, n), v(n, n), wrk(1)
      n2 = n + 1
      n3 = n + n2
      call svd1(a, m, n, m, n, 0, n, n, s, u, v, wrk, wrk(n2), wrk(n3))
      return 
      end
c  singular value decomposition of a matrix
c  Peter A. Businger, Bell Telephone Laboratories
c  Gene H. Golub, Stanford University
c  Algorithm 358, Collected algorithms from ACM
c  modified for real matrix by Mike Webring, USGS
      subroutine svd1(a, mmax, nmax, m, n, p, nu, nv, s, u, v, b, c, t)
      dimension a(mmax, 1), u(mmax, 1), v(nmax, 1)
      dimension s(n), b(n), c(n), t(n)
      integer p
      data eta / 1.5e-8 /
      data tol / 1.e-31 /
      np = n + p
c
c  householder reduction
      n1 = n + 1
      c(1) = 0.0
      k = 1
c
c  elimination of a(i,k), i=k+1...m
   10 k1 = k + 1
      z = 0.0
      do 20 i = k, m
   20 z = z + (a(i,k) ** 2)
      b(k) = 0.0
      if (z .le. tol) goto 70
      z = sqrt(z)
      b(k) = z
      w = abs(a(k,k))
      q = 1.0
      if (w .ne. 0.0) q = a(k,k) / w
      a(k,k) = q * (z + w)
      if (k .eq. np) goto 70
      do 50 j = k1, np
      q = 0.0
      do 30 i = k, m
   30 q = q + (a(i,k) * a(i,j))
      q = q / (z * (z + w))
      do 40 i = k, m
   40 a(i,j) = a(i,j) - (q * a(i,k))
c
c  phase transformation
   50 continue
      q = - sign(1.0,a(k,k))
      do 60 j = k1, np
c
c  elimination of a(k,j), j=k+2...n
   60 a(k,j) = q * a(k,j)
   70 if (k .eq. n) goto 140
      z = 0.0
      do 80 j = k1, n
   80 z = z + (a(k,j) ** 2)
      c(k1) = 0.0
      if (z .le. tol) goto 130
      z = sqrt(z)
      c(k1) = z
      w = abs(a(k,k1))
      q = 1.0
      if (w .ne. 0.0) q = a(k,k1) / w
      a(k,k1) = q * (z + w)
      do 110 i = k1, m
      q = 0.0
      do 90 j = k1, n
   90 q = q + (a(k,j) * a(i,j))
      q = q / (z * (z + w))
      do 100 j = k1, n
  100 a(i,j) = a(i,j) - (q * a(k,j))
c
c  phase transformation
  110 continue
      q = - sign(1.0,a(k,k1))
      do 120 i = k1, m
  120 a(i,k1) = a(i,k1) * q
  130 k = k1
c
c  tolerance for negligible elements
      goto 10
  140 eps = 0.0
      do 150 k = 1, n
      s(k) = b(k)
      t(k) = c(k)
  150 eps = amax1(eps,s(k) + t(k))
c
c  initialization of u and v
      eps = eps * eta
      if (nu .eq. 0) goto 180
      do 170 j = 1, nu
      do 160 i = 1, m
  160 u(i,j) = 0.0
  170 u(j,j) = 1.0
  180 if (nv .eq. 0) goto 210
      do 200 j = 1, nv
      do 190 i = 1, n
  190 v(i,j) = 0.0
c
c  qr diagonalization
  200 v(j,j) = 1.0
  210 do 380 kk = 1, n
c
c  test for split
      k = n1 - kk
  220 do 230 mm = 1, k
      m2 = (k + 1) - mm
      if (abs(t(m2)) .le. eps) goto 290
      if (abs(s(m2 - 1)) .le. eps) goto 240
c
c  cancellation of e(m2)
  230 continue
  240 cs = 0.0
      sn = 1.0
      m1 = m2 - 1
      do 280 i = m2, k
      f = sn * t(i)
      t(i) = cs * t(i)
      if (abs(f) .le. eps) goto 290
      h = s(i)
      w = sqrt((f * f) + (h * h))
      s(i) = w
      cs = h / w
      sn = - (f / w)
      if (nu .eq. 0) goto 260
      do 250 j = 1, n
      x = u(j,m1)
      y = u(j,i)
      u(j,m1) = (x * cs) + (y * sn)
  250 u(j,i) = (y * cs) - (x * sn)
  260 if (np .eq. n) goto 280
      do 270 j = n1, np
      q = a(m1,j)
      r = a(i,j)
      a(m1,j) = (q * cs) + (r * sn)
  270 a(i,j) = (r * cs) - (q * sn)
c
c  test for convergence
  280 continue
  290 w = s(k)
c
c  origin shift
      if (m2 .eq. k) goto 360
      x = s(m2)
      y = s(k - 1)
      g = t(k - 1)
      h = t(k)
      f = (((y - w) * (y + w)) + ((g - h) * (g + h))) / ((2.0 * h) * y)
      g = sqrt((f * f) + 1.0)
      if (f .lt. 0.0) g = - g
c
c  qr step
      f = (((x - w) * (x + w)) + (((y / (f + g)) - h) * h)) / x
      cs = 1.0
      sn = 1.0
      m1 = m2 + 1
      do 350 i = m1, k
      g = t(i)
      y = s(i)
      h = sn * g
      g = cs * g
      w = sqrt((h * h) + (f * f))
      t(i - 1) = w
      cs = f / w
      sn = h / w
      f = (x * cs) + (g * sn)
      g = (g * cs) - (x * sn)
      h = y * sn
      y = y * cs
      if (nv .eq. 0) goto 310
      do 300 j = 1, n
      x = v(j,i - 1)
      w = v(j,i)
      v(j,i - 1) = (x * cs) + (w * sn)
  300 v(j,i) = (w * cs) - (x * sn)
  310 w = sqrt((h * h) + (f * f))
      s(i - 1) = w
      cs = f / w
      sn = h / w
      f = (cs * g) + (sn * y)
      x = (cs * y) - (sn * g)
      if (nu .eq. 0) goto 330
      do 320 j = 1, n
      y = u(j,i - 1)
      w = u(j,i)
      u(j,i - 1) = (y * cs) + (w * sn)
  320 u(j,i) = (w * cs) - (y * sn)
  330 if (n .eq. np) goto 350
      do 340 j = n1, np
      q = a(i - 1,j)
      r = a(i,j)
      a(i - 1,j) = (q * cs) + (r * sn)
  340 a(i,j) = (r * cs) - (q * sn)
  350 continue
      t(m2) = 0.0
      t(k) = f
      s(k) = x
c
c  convergence
      goto 220
  360 if (w .ge. 0.0) goto 380
      s(k) = - w
      if (nv .eq. 0) goto 380
      do 370 j = 1, n
  370 v(j,k) = - v(j,k)
c
c  sort singular values
  380 continue
      do 450 k = 1, n
      g = -1.0
      j = k
      do 390 i = k, n
      if (s(i) .le. g) goto 390
      g = s(i)
      j = i
  390 continue
      if (j .eq. k) goto 450
      s(j) = s(k)
      s(k) = g
      if (nv .eq. 0) goto 410
      do 400 i = 1, n
      q = v(i,j)
      v(i,j) = v(i,k)
  400 v(i,k) = q
  410 if (nu .eq. 0) goto 430
      do 420 i = 1, n
      q = u(i,j)
      u(i,j) = u(i,k)
  420 u(i,k) = q
  430 if (n .eq. np) goto 450
      do 440 i = n1, np
      q = a(j,i)
      a(j,i) = a(k,i)
  440 a(k,i) = q
c
c  back transformation
  450 continue
      if (nu .eq. 0) goto 510
      do 500 kk = 1, n
      k = n1 - kk
      if (b(k) .eq. 0.0) goto 500
      q = - sign(1.0,a(k,k))
      do 460 j = 1, nu
  460 u(k,j) = q * u(k,j)
      do 490 j = 1, nu
      q = 0.0
      do 470 i = k, m
  470 q = q + (a(i,k) * u(i,j))
      q = q / (abs(a(k,k)) * b(k))
      do 480 i = k, m
  480 u(i,j) = u(i,j) - (q * a(i,k))
  490 continue
  500 continue
  510 if (nv .eq. 0) goto 570
      if (n .lt. 2) goto 570
      do 560 kk = 2, n
      k = n1 - kk
      k1 = k + 1
      if (c(k1) .eq. 0.0) goto 560
      q = - sign(1.0,a(k,k1))
      do 520 j = 1, nv
  520 v(k1,j) = q * v(k1,j)
      do 550 j = 1, nv
      q = 0.0
      do 530 i = k1, n
  530 q = q + (a(k,i) * v(i,j))
      q = q / (abs(a(k,k1)) * c(k1))
      do 540 i = k1, n
  540 v(i,j) = v(i,j) - (q * a(k,i))
  550 continue
  560 continue
  570 return 
      end
      function xty(x, ix, y, iy, n)
      dimension x(1), y(1)
      double precision tmp, tp, tn
      tp = 0.0d0
      tn = 0.0d0
      jx = 1
      jy = 1
      do 10 i = 1, n
      tmp = dble(x(jx)) * dble(y(jy))
      if (tmp .lt. 0.0d0) goto 5
      tp = tp + tmp
      goto 6
    5 tn = tn + tmp
    6 jx = jx + ix
   10 jy = jy + iy
      xty = sngl(tp + tn)
      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
      function ask_yn(request)
      character ask_yn*1
      character request*(*)
    2 write(unit=6, fmt=100) request
  100 format(1x,a,1h $)
      read(unit=5, fmt=102, err=1) ask_yn
  102 format(a1)
      if (ask_yn .eq. 'Y') ask_yn = 'y'
      if (ask_yn .eq. 'N') ask_yn = 'n'
      if ((ask_yn .ne. 'y') .and. (ask_yn .ne. 'n')) call scold(*2)
      return 
    1 call scold(*2)
c-----------------------------------------------------------------------
c--
      end
      subroutine scold(*)
      call bell
      write(unit=6, fmt=100) 
  100 format(/,38h *** Unacceptable response...Try again)
      return 1
c----------------------------------------------------------------------
c      subroutine ask(ques)
c      character*(*) ques
c      character*40 string
c      string=ques
c
c  --  This pads ques to left with blanks (or truncates).
c
c      write(6,101) string
c 101  format(5x,a40,$)
c      return
c      end
c-----------------------------------------------------------------------
c---
      end
      function correl(n, a, b)
      dimension a(n), b(n)
      correl = 0.0
      ab = 0.0
      a2 = 0.0
      b2 = 0.0
      do 1 i = 1, n
      ab = ab + (a(i) * b(i))
      a2 = a2 + (a(i) * a(i))
    1 b2 = b2 + (b(i) * b(i))
      a2b2 = a2 * b2
      if (a2b2 .lt. 1.e-30) return 
      correl = abs(ab) / sqrt(a2b2)
      return 
c***********************************************************************
c*******
      end
      function acosd2(dx, dh)
      acosd2 = 1.e+37
      hs = dh * dh
      xs = dx * dx
      if (hs .lt. xs) return 
      dy = sqrt(hs - xs)
      acosd2 = atan2(dy,dx) * 57.29578
      return 
c***********************************************************************
c*******
      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 if(rmax.ge.100.) then
        call xaxis(dxp, dyp, xp, 10., 2, .12, '(f5.0)', 4)
      else
        call xaxis(dxp, dyp, xp, 0.1, 2, .12, '(f4.1)', 4)
      endif
      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
c      subroutine home
c      character ding*1
c      ding = char(12)
c      write(unit=*, fmt=*) ding
c      return 
c      end
      subroutine sendesc(tail)
      character fmt*14, esc*1, tail*(*)
      n = len(tail)
      fmt = ' '
      assign 1 to label
      if (n .gt. 9) assign 2 to label
      if (n .gt. 99) assign 3 to label
      write(unit=fmt, fmt=label) n
    1 format(8h(1x,a1,a,i1,3h,$))
    2 format(8h(1x,a1,a,i2,3h,$))
    3 format(8h(1x,a1,a,i3,3h,$))
      esc = char(27)
      write(unit=6, fmt=fmt) esc, tail
      return 
      end


