重力勾配テンソル その2

続きです。

空間周波数fの定義は「回/m」。
対象範囲から最大・最小波長λを出して、最小・最大周波数fに変換すれば良いだけ。先の STACK OVERFLOW に python コードの例が載っていましたが、もっと簡潔に書けます。
周波数がわからなくても、波長から角波数k=2π/λを出せますし、波数だけならフーリエ変換時に０から順に出せるので、計算すら必要ないかもしれません。

曖昧な点は、実装しながら確認することに。

まずは文献に沿ったデータを用意。
次に「新・地震動のスペクトル解析入門」のソースを確認。これでも十分に短いのですが、python では numpy で1行、np.fft.fftn だけでOK。ありがたい。2D だと fft →転置→ fft →転置と同じ結果になったので、この順で内部処理しているだけかもしれません。

実装自体は容易で、2晩ほどでできました。

が、結果が合いません。似たような重力勾配や水平微分の分布になるのですが、オーダーの異なる場合があるのと、ky の符号が逆になっているように見えます。
角波数でダメならただの波数を試したり、ナイキスト周波数以上の領域を全て考慮したり、しなかったり。組み合わせによっては、それっぽい結果になるのですが、数値を記載している文献がないため確証を得られません。3日ほど考えましたが、最終的にはあきらめて寝かすことに。詳細に書かれた文献が、出てくるかもしれません。
もっと簡単に結果を出せると思っていたのですが、残念。

理解にはまだ時間が必要なようです。
ま、重力探査の現状と自分のレベルが分かっただけでも、良しとしましょう。

******************************************
20190813追記
fftn の並びに対応する fftfreq で使う符号、波数のあたりが怪しいようです。修正すると、それらしくなりました。が、検証できないのでここまで。もう少し寝かせましょう。