Gravmag prism fourier

[drandykass]

Here is the preliminary pull request for the prism Fourier calculations. This is wildly incomplete, but I wanted to put a pull request up so you could see the derivation and computations. Lots of messy comments right now.

I have a very peculiar issue that I'm working--the Fourier domain solution is offset from the space domain solution by (dx/2, dy/2). I THINK it has to do with how numpy computes the wavenumber vectors, but I'm not sure. Numpy returns a negative nyquist, which if not handled properly could result in this kind of shift in the phase spectrum. I thought everything was self consistent though. The following file demonstrates the issue.
utFourierGrav.pdf
If I write up everything manually, using my own wavenumber computation, folding routines, and only use the fft/ifft from Numpy, I don't think that error occurs.

Next steps (and thoughts):

1. Would it be prudent to have a sandbox directory that exists with the branch that contains testing files and such for collaboration, and then gets deleted from the branch before merging?
2. Figure out this dx/2 dy/2 issue.
3. Once this is complete, write up the rest of the gravity and magnetic functions (pretty low-hanging fruit).

You'll note that there is a general potential and a gravitational potential function. This is because both the mag and grav can use the general potential (1/r) formulation to construct their solutions.

Checklist:

• Make tests for new code
• Create/update docstrings
• Include relevant equations and citations in docstrings
• Code follows PEP8 style conventions
• Code and docs have been spellchecked
• Include new dependencies in docs, requirements.txt, README, and .travis.yml
• Documentation builds properly
• All tests pass
• Can be merged
• Changelog entry (leave for last)
• Firt-time contributor? Add yourself to doc/contributors.rst (leave for last)

[leouieda]

Thanks for openning this!

I have a very peculiar issue that I'm working--the Fourier domain solution is offset from the space domain solution by (dx/2, dy/2). I THINK it has to do with how numpy computes the wavenumber vectors, but I'm not sure. Numpy returns a negative nyquist, which if not handled properly could result in this kind of shift in the phase spectrum. I thought everything was self consistent though. The following file demonstrates the issue. utFourierGrav.pdf If I write up everything manually, using my own wavenumber computation, folding routines, and only use the fft/ifft from Numpy, I don't think that error occurs.

The FFT functions in every library are always a bit of a mistery. I have stumbled on these issues with the frequencies. What I came up with is in gravmag.transform. It would probably be better to make that function public (maybe even move it somewhere else?). Could you try using that to compute the frequencies and see if it works?

1. Would it be prudent to have a sandbox directory that exists with the branch that contains testing files and such for collaboration, and then gets deleted from the branch before merging?

Proably no need for that. The branch is itself a sandbox! You can put files anywhere you want and delete them after.

1. Once this is complete, write up the rest of the gravity and magnetic functions (pretty low-hanging fruit).

Oh, so you just make derivatives in the frequency domain of the potential? Cool! I never thought of that.

You'll note that there is a general potential and a gravitational potential function. This is because both the mag and grav can use the general potential (1/r) formulation to construct their solutions.

Ok, that's perfect. Maybe to keep with the convention, the gravitational potential could be just potential and the rest have different names?

One more thing, the way the gz function is now, it would not be a drop-in replacement for the space domain version. I think it would be better if the gz function receives input in the space domain and spits output also in the space domain, exactly the same as gravmag.prism. But we could have a gz_fft or gz_fd function that does the operation in the frequency domain. It would expect input in the frequency domain (kx, ky, kz, etc) and ouput in the frequency domain as well. Then, gz could do the FFT, pass it along to gz_fd and transform back again.

Having both functions would satisfy those wanting a drop-in replacement as well as those wanting to do more things in the frequency domain.

What do you think?