      subroutine adept(iw,icode,icode2,lx,xx,z,x,eps1,eps2)
c    estimates depth to magnetic basement from sampled magnetic profiles.
c    results must be interpreted using the convergence of the depth estimates.
c    subroutines required:  adept, depth, fork.
c    written by j.phillips.
c      common /magdata/ lx,xmin,xmax,amin,amax,zmin,zmax,
c     &xx(512),z(512),x(512)
      common /adep/a(5,512),r(5,512),v(512),d(4,512),nd(512),c
      common /work1/xc(512)
      common /work2/f(512)
      common /work3/b(512)
      common /work4/fb(512)
      dimension xx(lx),z(lx),x(lx)
      complex a,r,c,xc,f,b,fb,ftem,conjg,cm
      complex*16 xp
c      character*4 pri(2)/'p',' '/,sec(2)/'s',' '/
      character*1 sym
c      real x(512),d,title(18),aa(512),we(30),wo(29),fmt(18)
      real d,aa(512),we(30),wo(29)
      real*8 ap,xind
      data ncoef/16/
      data size/.08/
c      icode=100*(iw/2)
      delx=xx(2)-xx(1)
      xstart=xx(1)
c      print *,xstart,delx,lx
      nc=4
c      lxm1=lx-1
c     ap=-859.
      ap=0.
      do 5 i=1,lx
    5 ap=ap+x(i)
      ap=ap/lx
      do 11 i=1,lx
c      if(x(i).gt.amx) amx=x(i)
c      if(x(i).lt.amn) amn=x(i)
   11 xc(i)=x(i)-ap
c  determine prediction filter
      pi=3.1415927
      do 20 i=1,lx
      aa(i)=x(i)
   20 v(i)=x(i)
      we(1)=1.
      do 23 j=2,ncoef
      ap=0.
      xind=0.
      do 21 i=j,lx
      ap=ap+aa(i)*aa(i)+v(i-j+1)*v(i-j+1)
   21 xind=xind+aa(i)*v(i-j+1)
      rc=-2.*xind/ap
      do 22 i=j,lx
      temp=aa(i)
      aa(i)=aa(i)+rc*v(i-j+1)
   22 v(i-j+1)=v(i-j+1)+rc*temp
      we(j)=0.
      jh=(j+1)/2
      do 23 i=1,jh
      k=j-i+1
      temp=we(k)+rc*we(i)
      we(i)=we(i)+rc*we(k)
   23 we(k)=temp
c  extend profile
      k=lx+lx/2
      do 24 i=1,11
      lc=2**i
      if(lc.gt.k) go to 25
   24 continue
   25 if(lc.gt.512) then
        k=n
        do 26 i=1,11
        lc=2**i
        if(lc.gt.k) go to 27
   26   continue
   27   if(lc.gt.512) stop 'profile too long'
      endif
      lp=lx+1
      jh=(lp+lc)/2
      do 28 i=lp,jh
      xc(i)=0.
c      print *,lp,i,lc,lc-i+lp
      xc(lc-i+lp)=0.
      do 28 j=2,ncoef
      xc(i)=xc(i)-we(j)*xc(i-j+1)
   28 xc(lc-i+lp)=xc(lc-i+lp)-we(j)*xc(mod(lc-i+lp+j-2,lc)+1)
c  begin hilbert transform
      call fork(lc,xc,-1.)
      k=lc/2
      do 30 j=2,k
      xc(j)=xc(j)+xc(j)
   30 xc(lc+2-j)=0.
      call fork(lc,xc,+1.)
c  end hilbert transform
      do 40 i=1,lx
      v(i)=1.
      a(1,i)=1.
      r(1,i)=1.
      do 40 j=1,nc
   40 d(j,i)=0.
      is=iw/2
      iw=is+is+1
      iwm1=iw-1
      iwm2=iw-2
      iwm3=iw-3
      iwp1=iw+1
      isp1=is+1
      lxms=lx-is
      lmsp=lxms+1
c  initialize weights
      ze=4.
      do 41 j=1,is
      ro=j-is
      re=ro-0.5
      wo(j)=ro*ro+ze
      wo(iwm1-j)=wo(j)
      we(j)=re*re+ze
   41 we(iw-j)=we(j)
      dxs=ze
c  begin adaptive processing
      do 42 j=1,lx
      f(j)=xc(j)
   42 b(j)=xc(j)
c  enter main loop
      do 52 i=1,nc,2
      ip=i+1
      ib=(i-1)*is+2
      ie=lx-ib+i+1
      do 43 j=ib,ie
      fb(j)=f(j)*conjg(b(j-i))
   43 aa(j)=f(j)*conjg(f(j))+b(j-i)*conjg(b(j-i))
      ie=ie-iwm2
      do 48 k=ib,ie
      xp=0.
      ap=0.
      do 44 j=1,iwm2
      xp=xp+fb(k+j-1)/wo(j)
   44 ap=ap+aa(k+j-1)/wo(j)
      cm=(xp+xp)/ap
      xp=0.
      ap=0.
      do 45 j=1,iwm1
      xp=xp+fb(k+j-1)/we(j)
   45 ap=ap+aa(k+j-1)/we(j)
      c=(xp+xp)/ap
      ik=k+is-ip/2
      call depth(i,ik)
c  set weights
      zp=d(i,ik)*d(i,ik)
      if(zp.lt.dxs) zp=dxs
      do 46 j=1,iwm1
   46 we(j)=we(j)-ze+zp
      do 47 j=1,iwm2
   47 wo(j)=wo(j)-ze+zp
      ze=zp
      ik=ik+i/2
      ftem=f(ik)
      f(ik)=f(ik)-cm*b(ik-i)
   48 b(ik-i)=b(ik-i)-conjg(cm)*ftem
      ik=ik+ip/2
      f(ik)=f(ik)-c*b(ik-i)
      ib=ib+is
      ie=ie+is
      do 49 j=ib,ie
      fb(j)=f(j)*conjg(b(j-ip))
   49 aa(j)=f(j)*conjg(f(j))+b(j-ip)*conjg(b(j-ip))
      ie=ie-iwm3
      do 52 k=ib,ie
      xp=0.
      ap=0.
      do 50 j=1,iwm2
      xp=xp+fb(k+j-1)/wo(j)
   50 ap=ap+aa(k+j-1)/wo(j)
      c=(xp+xp)/ap
      ik=k+is-ip/2-1
      call depth(ip,ik)
c  set weights
      zp=d(ip,ik)*d(ip,ik)
      if(zp.lt.dxs) zp=dxs
      do 51 j=1,iwm2
   51 wo(j)=wo(j)-ze+zp
      ze=zp
      ik=ik+ip/2
      ftem=f(ik)
      f(ik)=f(ik)-c*b(ik-ip)
   52 b(ik-ip)=b(ik-ip)-conjg(c)*ftem
c  processing completed
c  begin output
c      sym='*'
      sym='+'
      do 60 i=1,is
      aa(i)=(i-1)*delx+xstart
   60 continue
c      d1=0.
c      d2=0.
c      v1=0.
c      v2=0.
      write(20,201) '   primary solutions'
      icon=0
      do 61 i=isp1,lxms
      aa(i)=(i-1)*delx+xstart
      j=nd(i)
      if(j.le.1) go to 61
      do 80 k=1,j
   80 d(k,i)=d(k,i)*delx
c      v(i)=3.*(d(1,i)-d(j,i))/(float(j-1)+.1e-05)
      v(i)=3.*(d(1,i)-d(j,i))/(float(j-1))
c      d0=d(1,i)
      d1=d(1,i)
c      v0=abs(v(i))
      v1=abs(v(i))
c      if(d0.gt.d1.and.d1.lt.d2) then
c      sym='.'
cc      if(v0.gt.v1.and.v1.lt.v2) then
          xm=aa(i)
          do 53 k=2,lx
          if(xm.gt.xx(k)) go to 53
          go to 54
   53     continue
   54     d1=d1+z(k-1)+(xm-xx(k-1))*(z(k)-z(k-1))/(xx(k)-xx(k-1))
c      if(v(i).gt.0.0) then
c        icode=400
c      else
c        icode=300
c      endif
      if(v1.lt.eps1) then
c        if(v1.gt.0.0) then
c        if(v(i).gt.0.0) then
cc          sym='+'
c          icode=600
c        else
cc          sym='-'
c          icode=200
c        endif
c      if(v1.lt.eps1) then
c      write(13) pri,aa(i-1),yy,d1-alt,v(i-1),
c     1 zero,zero,zero,zero
c          xm=aa(i-1)
c          call vchar(xm,d1,sym,1,icode,size,0.,0.,0.)
          write(20,202) xm,d1,v1
c       endif
c        if(icode.eq.300.or.icode.eq.400) then
c          icon=0
c          go to 61
c        endif
c          call line(xm,d1,1,0,icode)
          call line(xm,d1,1,icon,icode)
          icon=1
c          call line(xm,d1-v1,1,1,icode)
      else
        icon=0
        d1=z(k-1)+(xm-xx(k-1))*(z(k)-z(k-1))/(xx(k)-xx(k-1))
        if(v(i).gt.0.0) call vchar(xm,d1,sym,1,icode,size,0.,0.,0.)
      endif
c      d2=d1
c      d1=d0
c      v2=v1
c      v1=v0
   61 continue
  200 format(10x,3f7.3)
  201 format(a)
  202 format(5f15.6)
      do 62 i=lmsp,lx
      aa(i)=(i-1)*delx+xstart
   62 continue
c      print *,'62'
      its=1
c     write(11) its,(aa(i),xc(i),d(1,i),v(i),i=1,lx)
c  obtain implicit derivative result
c      sym='o'
c      sym='-'
      dxs=delx*delx
      isp1=is+is
      lxms=lx-isp1+1
c      d1=0.
c      d2=0.
c      v1=0.
c      v2=0.
      write(20,201) '   secondary solutions'
      icon=0
      do 75 k=isp1,lxms
      n=nd(k)
      if(n.le.1) go to 75
      do 70 i=1,n
   70 we(i)=d(i,k)*d(i,k)
c      print *,'70'
      d(n,k)=0.
      n=n-1
      do 71 i=1,n
      ip=i+1
      ip2=ip*ip
      i2=i*i
      arg=(ip2-i2+4.*(we(ip)-we(i))/dxs)
      if(arg.ne.0.0) arg=(ip2*we(i)-i2*we(ip))/arg
      if(arg.lt.0.) arg=0.
      d(i,k)=sqrt(arg)
   71 continue
c      print *,'71'
      aa(k)=(k-1)*delx+xstart
      iflag=1
      if(n.gt.2) go to 72
      v(k)=d(1,k)-d(n,k)
      go to 74
   72 v(k)=d(1,k)-d(2,k)
      nm1=n-1
      do 73 i=1,nm1
      ip1=i+1
      do 73 j=ip1,n
      cs=d(i,k)-d(j,k)
      if(abs(cs).gt.abs(v(k))) go to 73
      v(k)=cs
      iflag=i
   73 continue
   74 nd(k)=iflag
c      print *,k,nd(k),d(nd(k),k),v(k),'73'
c      d0=d(nd(k),k)
c      v0=abs(v(k))
      d1=d(nd(k),k)
          xm=aa(k)
          do 55 i=2,lx
          if(xm.gt.xx(i)) go to 55
          go to 56
   55     continue
   56     continue
c   56   print *,i,d1,xx(i),xx(i-1),'55'
c   56     d1=d1+z(i-1)+(xm-xx(i-1))*(z(i)-z(i-1))/(xx(i)-xx(i-1))
          d1=d1+z(i-1)+(xm-xx(i-1))*(z(i)-z(i-1))/(xx(i)-xx(i-1))
c      if(v(k).gt.0.0) then
c        icode=402
c      else
c        icode=302
c      endif
      v1=abs(v(k))
c      if(d0.gt.d1.and.d1.lt.d2) then
c      sym=','
cc      if(v0.gt.v1.and.v1.lt.v2) then
      if(v1.lt.eps2) then
c        if(v(k).gt.0.0) then
c          icode=602
c        else
c          icode=202
c        endif
c      print *,v1
c      write(14) sec,aa(k-1),yy,d1-alt,v(k-1),
c     1 zero,zero,zero,zero
c          xm=aa(k-1)
c          print *,xm,d1,icode,size
          write(20,202) xm,d1,v1
c      endif
c      if(icode.eq.302.or.icode.eq.402) then
c        icon=0
c        go to 75
c      endif
c          call line(xm,d1,1,0,icode)
c          call vchar(xm,d1,sym,1,icode,size,0.,0.,0.)
          call line(xm,d1,1,icon,icode2)
          icon=1
c          call line(xm,d1-v1,1,1,icode)
      else
        icon=0
        d1=z(i-1)+(xm-xx(i-1))*(z(i)-z(i-1))/(xx(i)-xx(i-1))
        if(v(i).gt.0.0) call vchar(xm,d1,sym,1,icode2,size,0.,0.,0.)
      endif
cc      endif
c      d2=d1
c      d1=d0
c      v2=v1
c      v1=v0
   75 continue
c      print *,'75'
      its=2
c     write(11) its,(aa(k),d(nd(k),k),v(k),k=isp1,lxms)
      its=-999
c     write(11) its
      return
      end
      subroutine depth(i,ik)
      common /adep/a(5,512),r(5,512),v(512),d(4,512),nd(512),c
      complex a,r,c,bot,conjg
c  determine new autocorrelation coefficient
      ip=i+1
      a(ip,ik)=0.
      r(ip,ik)=c*v(ik)
      v(ik)=v(ik)*(1.-c*conjg(c))
      do 20 j=2,ip
   20 r(ip,ik)=r(ip,ik)-a(j,ik)*r(ip-j+1,ik)
      ig=(ip+1)/2
      do 21 j=1,ig
      bot=a(ip-j+1,ik)-c*conjg(a(j,ik))
      a(j,ik)=a(j,ik)-c*conjg(a(ip-j+1,ik))
   21 a(ip-j+1,ik)=bot
c  determine depth
c     arg=(2*i-1)/(real(r(i,ik))/real(r(ip,ik))-1.)-(i-1)*(i-1)
      arg=real(r(ip,ik))
      if(arg.eq.0.0.or.arg.eq.1.0) go to 22
      arg=i*i/(1./arg-1.)
      if(arg.lt.0.) arg=0.
   22 d(i,ik)=0.5*sqrt(arg)
      nd(ik)=i
      return
      end

