feature: restyling of Kirchhoff operator

[mrava87]

This PR introduces a number of changes to the Kirchhoff operator in an attempt to improve the amplitude fidelity of this modelling process with respect to finite-difference full wavefield modelling. Also follows from #439.

Major changes:

• Remove anglerefl as the angle theta can be computed directly from the incident and reflection angles with normal
• Remove snell limitations, this is naturally handled by stationary-phase when summing along all possible image points and src-rec pairs
• Force amp to be provided as pair of tables like trav
• When amp is not provided, compute it as 1/sqrt(dist) for both source and receiver side in 2D (and 1/dist in 3D)
• Correct computation of cos(\theta)
• Documentation polished and fully consistent with code

[mrava87]

@cako, I got some time to revisit the Kirchhoff operator based on some the discussions we had in #439. The dynamic version looks much better now in terms of match with FD (see https://github.com/PyLops/pylops_notebooks/blob/master/developement/Modellingcomparison_kirch_vs_fd.ipynb).

Do you think you have some time to review this?

[cako]

I should have some time, yes. I will take a look!

[cako]

Great stuff! Code LGTM so I'm going to go ahead and approve it.

I do have a couple of comments/questions.

1. Why change the sign of the wavelet? Was it not matching the data as expected?
2. Applying the cosine of the reflection angle or using the apertures ensures some protection against very wide angles, but the Snell parameter is the best way to control this. However, upon rereading the removed code I don't think it was implemented properly so best to remove it for now. (E.g., if ang_src = 0 and ang_rec = 360, it would complain that the reflection angle is too big... it should be implemented on cos((ang_src - ang_rec)/2)).

[mrava87]

Great stuff! Code LGTM so I'm going to go ahead and approve it.

Thanks @cako, good you agree with the changes :)

I do have a couple of comments/questions.

1. Why change the sign of the wavelet? Was it not matching the data as expected?

If you refer to this Wfilt = W * (-1j * 2 * np.pi * f), I wasn't really sure myself but I think I made a mistake. If you look at eq 4 of Yang and Zhang, 2019 (which I added now as reference as I basically take all eqs from there) then it should indeed have a -. I fixed it now.

One thing I wasn't that sure is the 2D case, where in Safron they suggest using |\omega| but one of my visiting students who works on LSM suggested it should be \sqrt{jw} taken from Deregowski, S., and S. Brown, 1983, A theory of acoustic diffractors applied to 2d models: Geophysical Prospecting, 31, 293-333 and this seems to provide indeed better match with FD. Any comment?

2. Applying the cosine of the reflection angle or using the apertures ensures some protection against very wide angles, but the Snell parameter is the best way to control this. However, upon rereading the removed code I don't think it was implemented properly so best to remove it for now. (E.g., if ang_src = 0 and ang_rec = 360, it would complain that the reflection angle is too big... it should be implemented on cos((ang_src - ang_rec)/2)).

Indeed, I think that Snells code was wrong (actually all of my intuition over those angles was wrong). I see now why you just care about half-opening which is always implicitly suggesting reflection along he normal of the reflector even if you do not know what that is...

Regarding the cosine of the angle, I left it there as some formulations seems to suggest it should be there. However, what I noticed in Yang and Zhang, 2019 is that makes sense if you use r(x, theta) which is not what we effectively implement here (our r is only location dependent). Now if you make the simplification in their eq8 the angle gets absorbed. Not knowing where to lean, I decided to keep it as in some cases as you suggest you may want to use it as a weighting factor but if you do not want you can just force it like I do in the notebook (kop.angle_recs = kop.angle_srcs.copy() # done as a quick and dirty way to remove the cos(theta) weighting term in modelling). Again, any comment?

[cako]

Interesting! I didn't know that paper. Is there a significant difference between sqrt(jw) and |w|? In terms of frequency content they should be the same. The phase might be a bit altered. I think if yours matches FD simulation then it is probably the preferred way to go... with the caveat that traditionally Kirchhoff modelers try to do a good job at modeling the amplitude of reflections and not diffractions. A bit counter-intuitive in a way! In any case, I don't mind one way or the other, there are dozens of expressions in the literature, I think it's just a matter of choosing one.

Regarding the angle, I definitely agree with keeping the cos of the reflection angle. I just mean that one may want to further modify the diffraction stack based on a user-given weight. This can be useful for higher control when killing the tomographic term. In any case, we can handle it in a future PR if we want to tackle it.

Feel free to merge!