Necessary Solution: Horizontal loop source for r != 0

[dccowan]

We simulate the response due to a horizontal loop source by computing the solutions to Ward and Hohman equations 4.87 and 4.88

As of right now, we can only compute the solution for the radial distance r=0. In this case, the radial component of the fields is zero and the integral for the vertical component contains a single Bessel function (since J0(0) = 1). It would be great to have a stable solution that works everywhere.

Does empymod have the ability to design a filter that would compute these integrals in a stable manner?

[dccowan]

@sgkang @prisae I could use your input.

[prisae]

I have a couple of inputs here, but @sgkang might be of more help, as I don't know what (which formulation) is exactly implemented in simpegEM1D, and how general it is. Hopefully we can quickly touch on this in tomorrow's meeting as well.

1. The filter-design functionality (fdesign) within empymod is as such agnostic to the problem at hand, so you could definitely use it to design a filter if these two requirements are met:
   1. It is a linear transform (otherwise DLF won't work);
   2. You have an analytical transform pair, or a numerical transform pair with sufficiently accurate and precise results or a very wide range of wavenumbers.
2. You state "It would be great to have a stable solution that works everywhere." => However, then we might want to discuss if these equations are the right ones to implement. Be aware that these only hold if your large source loop is exactly horizontal, in parallel with the layers. If you want a general expression where your loop can also have an arbitrary orientation then you will need another formulation. What is often done is that you only compute the response to an infinitesimal small point dipole in the k-f domain, and take care of the source and receiver shape in the x-f domain (and of the waveform in the x-f or x-t domain). So to speak, you integrate a couple of point dipoles along your loop, usually with, e.g., Gauss-Legendre, it doesn't need many points. In the case where the loop actually is exactly horizontal it is also very cheap to do it in the frequency domain, as your field is symmetric around the z-axis, so you only need one kernel computation.
3. So it appears I am of no big help here... But I never used DLF for a kernel that had two Bessel functions, and I actually have never seen it being used in the literature. My take is that this equation is a very specific case, that has a nice integral representation for the x-f domain. However, numerically I would solve it differently.