subroutine newvec(grid,nx,ny,dx,dy,fi1,fd1,mi1,md1,fi2, & fd2,mi2,md2,store) c c Subroutine NEWVEC transforms a gridded total-field anomaly c into a new anomaly with new directions of magnetization and c ambient field. NEWVEC uses the following steps: (1) Fourier c transform the field, (2) multiply by the phase filter, and c (3) inverse Fourier transform the product. Anomaly values c are specified on a rectangular grid with x and y axes c directed north and east, respectively. Z axis is down. c Requires subroutines FOURN, DIRCOS, 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 total-field anomaly. 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 store - a singly dimensioned real array used internally. c It should be dimensioned at least 2*nx*ny in the c calling program. c dx - sample interval in the x direction, units irrelevant. c dy - sample interval in the y direction, units irrelevant. c mi1 - original inclination of magnetization, in degrees. c md1 - original declination of magnetization. c fi1 - original inclination of ambient field. c fd1 - original declination of ambient field. c mi2 - new inclination of magnetization, in degrees. c md2 - new declination of magnetization. c fi2 - new inclination of ambient field. c fd2 - new declination of ambient field. c c Output parameters: c grid - upon output, grid contains the transformed total- c field anomaly with same orientation as above. c dimension grid(nx*ny),store(2*nx*ny),nn(2) complex cgrid,cmplx,thetam1,thetam2,thetaf1,thetaf2, & cphase real kx,ky,k,mi1,md1,mi2,md2,mx1,my1,mz1,mx2,my2,mz2 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) call dircos(mi1,md1,0.,mx1,my1,mz1) call dircos(fi1,fd1,0.,fx1,fy1,fz1) call dircos(mi2,md2,0.,mx2,my2,mz2) call dircos(fi2,fd2,0.,fx2,fy2,fz2) do 10 j=1,nx do 10 i=1,ny ij=index(i,j,ny) store(2*ij-1)=grid(ij) 10 store(2*ij)=0. call fourn(store,nn,ndim,-1) do 20 j=1,nx do 20 i=1,ny ij=index(i,j,ny) if(ij.eq.1)then cphase=0. else call kvalue(i,j,nx,ny,dkx,dky,kx,ky) k=sqrt(kx**2+ky**2) thetam1=cmplx(mz1,(kx*mx1+ky*my1)/k) thetaf1=cmplx(fz1,(kx*fx1+ky*fy1)/k) thetam2=cmplx(mz2,(kx*mx2+ky*my2)/k) thetaf2=cmplx(fz2,(kx*fx2+ky*fy2)/k) cphase=thetam2*thetaf2/(thetam1*thetaf1) end if cgrid=cmplx(store(2*ij-1),store(2*ij)) cgrid=cgrid*cphase store(2*ij-1)=real(cgrid) 20 store(2*ij)=aimag(cgrid) call fourn(store,nn,ndim,+1) do 30 j=1,nx do 30 i=1,ny ij=index(i,j,ny) 30 grid(ij)=store(2*ij-1)/(nx*ny) return end