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