c GRADCOMP
c
c computes horizontal derivatives and horizontal gradient magnitude using
c the gradient-component method of Thurston and Brown (Geophysics, 1994).
c
c source code modified from Thurston and Brown (1994) by Jeff Phillips.
c
      parameter(maxnobs=21, ncmax=2000)
      real filtop1(maxnobs*maxnobs), filtop2(maxnobs*maxnobs)
      character ifile*80,id*56,pgm*8
      real surf(ncmax,maxnobs),grad(ncmax),xx(ncmax),yy(ncmax)
      data dval/1.e38/
c      open(12,file='filtop1.xyz',form='unformatted',status='unknown')
c      open(13,file='filtop2.xyz',form='unformatted',status='unknown')
c get the grid to be filtered
      print*,'enter grid to be filtered'
      read(5,100) ifile
  100 format(a)
      open(10,file=ifile,status='old',form='unformatted')
      print*,'enter output horizontal gradient grid'
      read(5,100) ifile
      open(11,file=ifile,status='unknown',form='unformatted')
      print*,'enter output x-derivative grid'
      read(5,100) ifile
      open(12,file=ifile,status='unknown',form='unformatted')
      print*,'enter output y-derivative grid'
      read(5,100) ifile
      open(13,file=ifile,status='unknown',form='unformatted')
      read(10) id,pgm,nc,nr,nz,xo,dx,yo,dy
      write(11) id,pgm,nc,nr,nz,xo,dx,yo,dy
      write(12) id,pgm,nc,nr,nz,xo,dx,yo,dy
      write(13) id,pgm,nc,nr,nz,xo,dx,yo,dy
c input number of points in operator and polynomial order
   31 print*,'enter the # if weights in the filter (odd and < 22)'
      read(5,*) nobs
      odd=mod(nobs,2)
      if(odd.eq.0) then
        print*,'re-enter as odd'
        go to 31
      endif
      if(nobs.gt.21) then
        print*,'must be less than 22'
        go to 31
      endif
   32 print*,'input order of polynomial'
      read(5,*) order
      if(order.ge.nobs) then
        print*,'order must be less than window size'
        go to 32
      endif
      call gradop(filtop1,filtop2,nobs,order,dx)
c      do 10 i=1,nobs
c      y=i
c      do 10 j=1,nobs
c      x=j
c      write(12) x,y,filtop1(j+(i-1)*nobs)
c   10 continue
c      do 11 i=1,nobs
c      y=i
c      do 11 j=1,nobs
c      x=j
c      write(13) x,y,filtop2(j+(i-1)*nobs)
c   11 continue
      mid=(nobs-1)/2+1
      do i=1,nc
      grad(i)=dval
      xx(i)=dval
      yy(i)=dval
      enddo
      do j=1,mid-1
      write(11) dlt,(grad(i),i=1,nc)
      write(12) dlt,(xx(i),i=1,nc)
      write(13) dlt,(yy(i),i=1,nc)
      enddo
      do j=1,nobs-1
      read(10) dlt,(surf(i,j),i=1,nc)
      enddo
      do k=mid,nr-mid+1
      read(10) dlt,(surf(i,nobs),i=1,nc)
      call convolve(surf,grad,xx,yy,filtop1,filtop2,nc,ncmax,nobs)
      write(11) dlt,(grad(i),i=1,nc)
      write(12) dlt,(xx(i),i=1,nc)
      write(13) dlt,(yy(i),i=1,nc)
      do j=1,nobs-1
      do i=1,nc
      surf(i,j)=surf(i,j+1)
      enddo
      enddo
      enddo
      do i=1,nc
      grad(i)=dval
      enddo
      do j=nr-mid+2,nr
      write(11) dlt,(grad(i),i=1,nc)
      write(12) dlt,(grad(i),i=1,nc)
      write(13) dlt,(grad(i),i=1,nc)
      enddo
      stop
      end
c
      subroutine gradop(filtop1,filtop2,nobs,order,grdint)
c
c This subroutine generates two gradient-component operators
c using the method from Thurston and Brown (Geophysics, 1994)
c Input variable:
c    nobs   - size of gradient-component operator
c Output variables:
c    filtop1 - d/dx gradient-component operator
c    filtop2 - d/dy gradient component operator
c Internal variables:
c    orpol   - array of 1d orthogonal polynomials stored as
c              a 2d array according to order and abscissa
c    dorpol  - d/dx of orpol
c    den     - sum of squares of orthogonal polynomials
c
      parameter (maxord=20,maxnobs=21)
      real orpol(-2:maxord,maxnobs),dorpol(-2:maxord,maxnobs)
      real den(0:maxord), filt(maxnobs), dfilt(maxnobs)
      real filtop1(nobs,nobs), filtop2(nobs,nobs)
      mid=(nobs+1)/2
c calculate orthogonal polynomials and derivatives
      call orpoly(order,nobs,orpol,dorpol)
c calculate gradient-component operators
      do i=0,order
      do j=1,nobs
        den(i)=0.
      enddo
      enddo
c calculate denominator
      do i=0,order
      do j=1,nobs
        den(i)=den(i)+orpol(i,j)**2
      enddo
      enddo
c calculate 1-D smoothing and derivative operators
      do 10 i=1,nobs
        filt(i)=0.
        dfilt(i)=0.
      do 11 j=0,order
        filt(i)=filt(i)+orpol(j,i)*orpol(j,mid)/den(j)
        dfilt(i)=dfilt(i)+orpol(j,i)*dorpol(j,mid)/den(j)
   11 continue
   10 continue
c combine 1-D operators to generate 2-D operators
      do 15 i=1,nobs
      do 16 j=1,nobs
        filtop1(i,j)=filt(j)*dfilt(i)/grdint
        filtop2(i,j)=filt(i)*dfilt(j)/grdint
   16 continue
   15 continue
      return
      end
c
      subroutine orpoly(order,nobs,orpol,dorpol)
c this subroutine generates orthogonal polynomials in x, and
c their first derivatives wrt x.
      parameter(maxord=20,maxnobs=21)
      real orpol(-2:maxord,maxnobs), dorpol(-2:maxord,maxnobs)
      real ssqpol(0:maxord),x(maxnobs)
      integer t
c initialize array containing x
      do 8 i=1,nobs
        x(i)=i
    8 continue
c initialize the first two elements of orpol and dorpol for all x
      do 1 i=1,nobs
        orpol(-2,i)=-1.
        orpol(-1,i)=0.
        dorpol(-2,i)=0.
        dorpol(-1,i)=0.
    1 continue
c generate successive polynomials using recursion
      alpha=0.0
      beta=1.0
      oldsqsum=1.0
      do 2 t=0,order
        ssqpol(t)=0.
        xprod=0.
        do 3 i=1,nobs
          orpol(t,i)=orpol(t-2,i)*beta + orpol(t-1,i)*(x(i)+alpha)
          dorpol(t,i)=dorpol(t-2,i)*beta + dorpol(t-1,i)*(x(i)+alpha)
     +               + orpol(t-1,i)
          d=orpol(t,i)**2
          ssqpol(t)=ssqpol(t)+d
          xprod=xprod+d*x(i)
    3   continue
        beta=-ssqpol(t)/oldsqsum
        oldsqsum=ssqpol(t)
        alpha=-xprod/oldsqsum
    2 continue
      return
      end
c
      subroutine convolve(surf,grad,xx,yy,filtop1,filtop2,nc,ncmax,nobs)
      real surf(ncmax,nobs),grad(nc),xx(nc),yy(nc)
      real filtop1(nobs,nobs), filtop2(nobs,nobs)
      data ddval/1.e30/,dval/1.e38/
      mid=(nobs-1)/2+1
c      do 10 k=1,mid-1
c      grad(k)=dval
c      grad(nc-k+1)=dval
c   10 continue
      do 70 k=mid,nc-mid+1
        do 40 j=1,nobs
        do 40 i=1,nobs
          if(surf(i+k-mid,j).ge.ddval) then
            grad(k)=dval
            xx(k)=dval
            yy(k)=dval
            go to 70
          endif
   40   continue
        grad1=0.
        grad2=0.
        do 50 j=1,nobs
        do 50 i=1,nobs
          grad1=grad1+surf(i+k-mid,j)*filtop1(i,j)
          grad2=grad2+surf(i+k-mid,j)*filtop2(i,j)
   50   continue
        grad(k)=sqrt(grad1**2+grad2**2)
        xx(k)=grad1
        yy(k)=grad2
   70   continue
        return
        end
