Add potential/gravprisms as new module for computing geopotential anomalies

[PaulWessel]

Description of proposed changes

This PR adds gravprisms which complements recent improvements to grdseamount related to a variable density structure:

1. It can compute the faa, geoid, or vgg over a set of vertical prisms - these may have different density values
2. Instead of reading prisms it can create prisms from a seamount density model given via grids.

The documentation to gravprisms includes a new figure created with it and plot3d, I show it here:

I have not added any test scripts yet (will try to add numerical scripts) but I have compared its output to other GMT modules for the cases we can compare and it is solid:

For constant density between two surfaces it compares well with the more exact output from grdgravmag3d but is much faster. For variable rho(r) is is more accurate than gravfft, and for rho(r,z) it is our only option.

I have not yet tested the vgg and geoid output hence WIP.

[PaulWessel]

I have now added a single numerical tests that checks that the output for faa, vgg, and geoid is correct.

[PaulWessel]

This module is now ready for review. The single numerical test passes for gravity, VGG, and geoid and the output has been vetted with other modules. I have added checks for the somewhat unstable geoid expressions (48 terms per prism) to catch exceptions (where a component should be zero rather than calling atan and log with 1/0 and 0 as arguments, respectively). Documentation mentions this regarding the geoid. I have done lots of comparisons with gravfft, talwani3d and grdgravmag3d for the cases that can be compared and all looks very good - exact matches for constant density cases.

[joa-quim]

Sorry for bringing this up so late but maybe rename this module to gravmagprism. I know It doesn't (yet) do mag case but it should only be a matter of replacing the grav prism function for the mag one, which can be done later.

[PaulWessel]

Not too late at all. The only issue (not specific to this one) is "gravmag". I guess we are stuck with that but geopotential would have been a better thing even if long....

So gravmagprism it is. I can add language that says magnetics not yet implemented.

Please tell me what options we will need in the future to do magnetics. SOmething like -H from grdgravmag3d.c?

These options are available still: -Q, as well as -K, -P since not a plotter. I am using quite a few options to handle the many inputs so that it would be simple to use from externals...

[PaulWessel]

ANother thought: We still have just gravfft even though one can do something very similar for magnetics, no? Given the lack of options available and the possibility of needing specific magnetic options, perhaps gravprisms is not so bad - and when you get energized you can duplicate it to build magprisms.c?

[joa-quim]

Yeap, a eventual future magprisms is OK. For the moment grdgravmag3d is resisting like hell to give a convincing result with variable density.

[PaulWessel]

For the moment grdgravmag3d is resisting like hell to give a convincing result with variable density.

You have to be strict and show it who's the boss.

[joa-quim]

There is a previous problem in the two layers case and that's why it fails, I think.

[maxrjones]

@PaulWessel do you have the prisms.txt file used in the examples for this module? If so, could we make this a remote file so that the examples are reproducible?

[PaulWessel]

You mean for the plots at the top of this issue? Or the stuff under EXAMPLES in the docs?

[maxrjones]

The stuff under EXAMPLES in the docs.