      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
