subroutine signal(grid,nx,ny,dx,dy,store1,store2,store3)
c
c  Subroutine SIGNAL  calculates the analytical signal of 
c  gridded potential field data using the following steps:  
c  (1) Fourier transform the field and make two copies, 
c  (2) multiply these arrays by ikx, iky, and k, respectively, 
c  (3) inverse Fourier transform all three arrays, and (4) 
c  calculate the square root of the sum of the squares of the 
c  three arrays.   Field values are specified on a rectangular grid 
c  with x, y, and z axes directed north, east, and down, 
c  respectively;  North is arbitrary.  Requires subroutines 
c  FOURN and KVALUE. 
c
c  Input parameters:
c    nx - number of elements in the south-to-north direction.  
c    ny - number of elements in the west-to-east direction.  
c         (NOTE:  both nx and ny must be a power of two.)
c    grid - a singly dimensioned real array containing the 
c         two-dimensional potential field.  Elements should 
c         be in order of west to east, then south to north (i.e., 
c         element 1 is southwest corner, element ny is 
c         southeast corner, element (nx-1)*ny+1 is northwest 
c         corner, and element ny*nx is northeast corner. 
c    store1,store2,store3 - singly dimensioned real arrays used 
c         internally.  Each should be dimensioned at least 2*nx*ny
c         in the calling program. 
c    dx - sample interval in the x direction, units irrelevant.
c    dy - sample interval in the y direction, units irrelevant.
c
c  Output parameters:
c    grid - upon output, grid contains the analystical signal of
c         the potential field with same orientation as above.
c
      dimension grid(nx*ny),store1(2*nx*ny),store2(2*nx*ny),
     &          store3(2*nx*ny),nn(2)
      complex cgrid1,cgrid2,cgrid3,cmplx
      real kx,ky,k
      data pi/3.14159265/
      index(i,j,ncol)=(j-1)*ncol+i
      nn(1)=ny
      nn(2)=nx
      ndim=2
      dkx=2.*pi/(nx*dx)
      dky=2.*pi/(ny*dy)
      do 10 j=1,nx
         do 10 i=1,ny
            ij=index(i,j,ny)
            store1(2*ij-1)=grid(ij)
   10       store1(2*ij)=0.
      call fourn(store1,nn,ndim,-1)
      do 40 i=1,2*nx*ny
         store2(i)=store1(i)
   40    store3(i)=store1(i)
       do 20 j=1,nx
         do 20 i=1,ny
            ij=index(i,j,ny)
            call kvalue(i,j,nx,ny,dkx,dky,kx,ky)
            k=sqrt(kx**2+ky**2)
            cgrid1=cmplx(store1(2*ij-1),store1(2*ij))
            cgrid2=cmplx(store2(2*ij-1),store2(2*ij))
            cgrid3=cmplx(store3(2*ij-1),store3(2*ij))
            cgrid1=cgrid1*cmplx(0,kx)
            cgrid2=cgrid2*cmplx(0,ky)
            cgrid3=cgrid3*k
            store1(2*ij-1)=real(cgrid1)
            store1(2*ij)=aimag(cgrid1)
            store2(2*ij-1)=real(cgrid2)
            store2(2*ij)=aimag(cgrid2)
            store3(2*ij-1)=real(cgrid3)
            store3(2*ij)=aimag(cgrid3)
   20       continue
      call fourn(store1,nn,ndim,+1)
      call fourn(store2,nn,ndim,+1)
      call fourn(store3,nn,ndim,+1)
      do 30 j=1,nx
         do 30 i=1,ny
            ij=index(i,j,ny)
            grid1=store1(2*ij-1)/(nx*ny)
            grid2=store2(2*ij-1)/(nx*ny)
            grid3=store3(2*ij-1)/(nx*ny)
   30       grid(ij)=sqrt(grid1**2+grid2**2+grid3**2)
      return
      end