      subroutine hilbert(iopt,ipen,sym,xmin,xmax,zmin,zmax,
     1ainc,adec,n,a,x,z)
c
c  Depth, dip, and magnetic susceptibility estimates from total-field magnetic
c  profile data using the analytic signal method of Nabighian (1972, 1974)
c  (Geophysics, v.37 no.3, p.507-517; and Geophysics, v.39 no.1, p.85-92).
c
c  On input a(n) contains either the magnetic profile (for sheet solutions) or
c  the horizontal derivative of the magnetic profile (for contact solutions).
c  x(n) and z(n) contain the horizontal and vertical coordinates of the
c  observations.  a(n), x(n), and z(n) must be evenly sampled and distance
c  units must be consistent.  Magnetic inclination (ainc) and declination
c  (adec) are required to estimate dip and magnetic susceptibility.  The pen
c  number (ipen), symbol number (sym), x-range (xmin, xmax) and depth range
c  (zmin, zmax) are used for plotting.  To plot the analytic signal amplitude
c  curves, set iopt = 2, otherwise set iopt = 0.
c
      common /work1/xc(512)
      common /work2/hx(512),vx(512)
      common /work3/rpx(512),tmp(512)
      common /model/nbod,ncor,iopn,delx,ibod
      common /magxyz/ ef(3),ev(3),tdir(3),azmuth,idum
      dimension a(n),x(n),z(n),dip(3),susc(3)
      dimension aa(3),bb(3)
      complex xc
      character sym*1
        data d2r/1.745329e-2/
c
c  Compute the complex analytic signal xc(lc) of the input profile a(n)
c
      size=.08
      call setfft(1,lc,n,a)
      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
c  Compute (and plot) the analytic signal amplitude tmp(n)
c
      tmax=-1.e+38
      do 40 j=1,n
      rpx(j)=0.
      tmp(j)=real(xc(j))**2+aimag(xc(j))**2
      if(x(j).ge.xmin.and.x(j).le.xmax) then
        if(tmp(j).gt.tmax) tmax=tmp(j)
      endif
   40 continue
      tmin=0.0
      if(iopt.gt.1) then
        call dplot(2,xmin,xmax,tmin,tmax)
        call curv(x,tmp,n,ipen+2)
      endif
      call dplot(1,xmin,xmax,zmin,zmax)
c
c  Set up constants
c
      dx=x(2)-x(1)
      if(ainc.lt.1.e38) then
        capa=(azmuth-adec)*d2r
        term=2.*atan2(tan(ainc*d2r),cos(capa))/d2r
        clc=-2.*(1.-(cos(ainc*d2r)*sin(capa))**2)
      endif
c
c  Loop to locate maxima along the analytic signal amplitude curve
c
      do 45 j=2,n-1
      call interp(tmp(j-1),tmp(j),tmp(j+1),dx,xo,yo,ierr)
      if(ierr.ne.0) go to 45
c
c  Maximum located, get 3 estimates of (squared) depth and use the minimum one
c
      z1=tmp(j-1)*(-dx-xo)**2/(yo-tmp(j-1))
      if(yo.ne.tmp(j)) then
        z2=tmp(j)*(-xo)**2/(yo-tmp(j))
      else
        z2=z1
      endif
      z3=tmp(j+1)*(dx-xo)**2/(yo-tmp(j+1))
      z1=amin1(z1,z2,z3)
      xo=x(j)+xo
   41 continue
      z1=sqrt(z1)
      z1sav=z1
c
c  Correct the depth for the observation altitude and plot it
c
      z2=z1+z(j-1)+(xo-x(j-1))*(z(j+1)-z(j-1))/(x(j+1)-x(j-1))
      call vchar(xo,z2,sym,1,ipen,size,0.,0.,0.)
      if(ainc.ge.1.e38) go to 45
c
c  Get 3 estimates for the dip and magnetic susceptibility
c
      ak=z1sav*sqrt(yo)
      do 11 i=j-1,j+1
      k=i-j
      kk=k+2
      fz=real(xc(i))
      fx=aimag(xc(i))
      phi=((x(i)-xo)*fx-z1sav*fz)/(z1sav*fx+(x(i)-xo)*fz)
      phi=atan(phi)/d2r
      dip(kk)=term-phi
        if(dip(kk).lt.0.) then
          dip(kk)=dip(kk)+180.
        endif
        if(dip(kk).gt.180.) then
          dip(kk)=dip(kk)-180.
        endif
          dr=dip(kk)*d2r
       susc(kk)=-sign(ak,fz/cos(d2r*phi))/(clc*ef(1)*sin(dr))
c
c  Plot the three dip estimates
c
        capa=(zmax-zmin)/20.
            aa(1)=xo
            aa(2)=aa(1)+capa*cos(dr)
            bb(1)=z2
            bb(2)=bb(1)+capa*sin(dr)
            call curv(aa,bb,2,ipen)
11    continue
c
c  Write the location, depth, dips, and susceptibilities to 'pdepth.dep'
c
            write(20,201) xo,z2,dip(1),dip(2),dip(3),
     1susc(1),susc(2),susc(3)
201         format(5f15.6,/,30x,3f15.6)
   45 continue
   46 if(iopt.ne.1.and.iopt.ne.3) return
c
c  The rest of this is no longer used
c
      size=.06
      do 50 j=1,n
      vx(j)=real(xc(j))
      rpx(j)=aimag(xc(j))
   50 continue
      call spline(vx,hx,n,1.)
      call spline(rpx,vx,n,1.)
      tmin=1.e+38
      tmax=-1.e+38
      do 55 j=1,n
      tmp(j)=(hx(j)**2+vx(j)**2)
      if(x(j).ge.xmin.and.x(j).le.xmax) then
        if(tmp(j).lt.tmin) tmin=tmp(j)
        if(tmp(j).gt.tmax) tmax=tmp(j)
      endif
   55 continue
      if(iopt.gt.1) then
        call dplot(2,xmin,xmax,tmin,tmax)
        call curv(x,tmp,n,ipen+3)
      endif
      call dplot(1,xmin,xmax,zmin,zmax)
      do 60 j=2,n-1
      call interp(tmp(j-1),tmp(j),tmp(j+1),dx,xo,yo,ierr)
      if(ierr.ne.0) go to 60
      z3=1.e+38
      aa=tmp(j-1)-yo
      b=2.*tmp(j-1)*(-dx-xo)**2
      c=tmp(j-1)*(-dx-xo)**4
      call quadrat(aa,b,c,z1,z2,ierr)
      go to (63,62,61) ierr+1
   61 z3=amin1(z1,z2)
      go to 64
   62 z3=z1
   63 aa=tmp(j)-yo
      b=2.*tmp(j)*(xo)**2
      c=tmp(j)*(xo**4)
      call quadrat(aa,b,c,z1,z2,ierr)
      go to (66,65,64) ierr+1
   64 z3=amin1(z1,z2,z3)
      go to 64
   65 z3=amin1(z1,z3)
   66 aa=tmp(j+1)-yo
      b=2.*tmp(j+1)*(dx-xo)**2
      c=tmp(j+1)*(dx-xo)**4
      call quadrat(aa,b,c,z1,z2,ierr)
      go to (69,68,67) ierr+1
   67 z3=amin1(z1,z2,z3)
      go to 64
   68 z3=amin1(z1,z3)
   69 if(z3.eq.1.e+38) go to 60
      z3=sqrt(z3)
      xo=x(j)+xo
      if(abs(x(j+1)-x(j-1)).lt.1.e-30) go to 60
      z3=z3+z(j-1)+(xo-x(j-1))*(z(j+1)-z(j-1))/(x(j+1)-x(j-1))
      call vchar(xo,z3,sym,1,ipen,size,0.,0.,0.)
   60 continue
      return
      end
c---------------------------------------------------------------------------
      subroutine interp(y1,y2,y3,d,xmax,ymax,ierr)
c
c     Subroutine interp fits a parabola through ordinants y1,y2 and y3, assumed
c  to be evenly spaced, and reports the x and y coordinate of the parabola's
c  maximum.  Note...y1 must be .lt. y2 and y3 must be .lt. y2.
c
      ierr=1
      if(y1.lt.y2.and.y3.lt.y2) then
        ierr=0
      else
        return
      endif
      if(abs(d).lt.1.e-30) stop 'd'
      a=0.5*(y1-2.*y2+y3)/(d*d)
      b=0.5*(y3-y1)/d
      c=y2
      if(a.eq.0.)then
        ierr=1
        return
      end if
      if(abs(a).lt.1.e-30) stop 'a'
      xmax=-0.5*b/a
      ymax=a*xmax**2+b*xmax+c
      return
      end
c******************************************************************************
      subroutine quadrat(a,b,c,z1,z2,isol)
      isol=0
      d=b**2-4.*a*c
      if(d.lt.0.) go to 10
      if(abs(a).lt.1.e-30) go to 10
      d=sqrt(d)
      a=2.*a
      isol=2
      z1=(-b+d)/a
      z2=(-b-d)/a
      if(z2.lt.0.) isol=isol-1
      if(z1.lt.0.) then
        z1=z2
        isol=isol-1
      endif
   10 return
      end
c******************************************************************************
      subroutine eval(a,b,c,value)
      dimension a(1),b(1),c(1)
      value=a(1)*(b(2)*c(3)-b(3)*c(2))-b(1)*(a(2)*c(3)-a(3)*c(2))
     & +c(1)*(a(2)*b(3)-a(3)*b(2))
      return
      end
