[REVIEW] Extend documentation for power spectrum methods

[santisoler]

Part of https://github.com/openjournals/joss-reviews/issues/1544.

CuriedGrid has the following methods: radial_spectrum and radial_spectrum_log. Their names suggest the second one computes the radial spectrum and applies the np.log, while the first one does not. Nevertheless, the difference between both of them does not rely on the np.log but on the np.sqrt:

CurieGrid.radial_spectrum https://github.com/brmather/pycurious/blob/759fc520c22503f0b40b6a63faa4e632dd9f1e88/pycurious/grid.py#L250

CurieGrid.radial_spectrum_log https://github.com/brmather/pycurious/blob/759fc520c22503f0b40b6a63faa4e632dd9f1e88/pycurious/grid.py#L318

Moreover, I don't quite understand what's the goal of CurieGrid.azimuthal_spectrum. That could be related to some lack of knowledge on the subject from my part.

I think extending their docstring would help to a better understanding of what they are doing. Maybe including some math and references would be nice. Renaming the methods in order to make very clear what they do it's also a good approach.

Am I understanding something wrong here or do you agree with my point of view?

[brmather]

Thank you for raising this issue @santisoler - it is good to revisit these CurieGrid methods with fresh :eyes: because radial_spectrum and radial_spectrum_log are, as you say, identical except for one line. The difference comes down to slight differences between the Bouligand et al. (2009) and Tanaka et al. (1999) methods to compute Curie depth - i.e. squaring the FFT of the magnetic anomaly versus taking the square root.

It is helpful to understand that

2*log(FT) == log(FT**2)

thus I was able to refactor the code in #19 to pass a constant parameter to control whether we are taking the square or square root:

https://github.com/brmather/pycurious/blob/cff769036d73905629a448812853614afda611b3/pycurious/grid.py#L269

Setting const=0.5 is the same as np.log(np.sqrt(FT[mask])).

[brmather]

Regarding azimuthal_spectrum, this method is essentially used to polarise the magnetic anomaly for a range of azimuths. The azimuths are binned according to the angle increment theta (in degrees). So you essentially retrieve a stack of power spectra that correspond to different angles of polarisation.

@rdelhaye perhaps you would like to elaborate on this and improve the docstring for azimuthal_spectrum and also at the end of notebook 1. While we didn't get very far, if I remember correctly you wanted to use this to estimate Curie depth in regions that exhibit a strong magnetic orientation...? In any case we should put more maths in our docstrings.

[santisoler]

Thank you for raising this issue @santisoler - it is good to revisit these CurieGrid methods with fresh eyes because radial_spectrum and radial_spectrum_log are, as you say, identical except for one line. The difference comes down to slight differences between the Bouligand et al. (2009) and Tanaka et al. (1999) methods to compute Curie depth - i.e. squaring the FFT of the magnetic anomaly versus taking the square root.

It is helpful to understand that

2*log(FT) == log(FT**2)

thus I was able to refactor the code in #19 to pass a constant parameter to control whether we are taking the square or square root:

https://github.com/brmather/pycurious/blob/cff769036d73905629a448812853614afda611b3/pycurious/grid.py#L269

Setting const=0.5 is the same as np.log(np.sqrt(FT[mask])).

My mind was on "coder" mode while reading this and haven't figured out that the main difference is the power of the FT inside the np.log. Thanks for bringing my mind back to "math" mode haha.

In fact, your proposal of setting a const is a great choice!

[rdelhaye]

@brmather Thanks for looping me in on this! The intent of the radial_spectrum function, as Ben mentioned, was to examine the isotropy of the magnetic data.

One of the assumptions commonly made in literature with the Tanaka-style spectral analysis is that the magnetisation is essentially isotropic, with the spatial Fourier transforms averaged with respect to azimuth. We thought that having the ability to instead examine the transforms on a more granular basis, by subdividing the transforms into bins by azimuth, would lead to at least the opportunity to test how valid the assumption is. I'm not entirely certain if spectral analysis results in an anisotropic area are a) valid and b) useful (i.e., which value to use? The radially averaged results? Or the deeper results? I would say that the shallower result, based on the assumptions that we're computing the bottom, is probably "more wrong").

[brmather]

I think the key point to make here is that we leave open the possibility for the user to examine anisotropy of the magnetic field. In the future we can look at implementing additional Jupyter notebooks to explore this in greater detail, but for version 1.0 we have all the ingredients to compute Curie depth using two different approaches.

[brmather]

CuriedGrid has the following methods: radial_spectrum and radial_spectrum_log. Their names suggest the second one computes the radial spectrum and applies the np.log, while the first one does not. Nevertheless, the difference between both of them does not rely on the np.log but on the np.sqrt

I have submitted pull request #26 to eliminate radial_spectrum_log from the CurieGrid class and use radial_spectrum with the power argument to toggle between Bouligand or Tanaka approaches to compute Curie depth. Docstrings have been revised alongside Jupyter notebooks. Feel free to review the changes before I merge into master.