.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples/plot_internal_external.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_plot_internal_external.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_plot_internal_external.py:


Internal / external field separation
------------------------------------

Ground magnetometers see both the ionospheric (external) currents and
the telluric currents induced inside the Earth (internal). Placing
divergence-free SECS on two shells -- one in the ionosphere and one
below the ground -- separates the two contributions, as suggested in
the conclusions of Amm and Viljanen (1999). All three magnetic field
components are needed for the separation.

.. GENERATED FROM PYTHON SOURCE LINES 12-78



.. image-sg:: /examples/images/sphx_glr_plot_internal_external_001.png
   :alt: Two-shell separation of external and internal sources
   :srcset: /examples/images/sphx_glr_plot_internal_external_001.png
   :class: sphx-glr-single-img





.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np

    from pysecs import SECS, T_df


    R_E = 6378e3
    rng = np.random.default_rng(5)

    # External shell in the ionosphere, internal image shell below ground
    lat_g, lon_g = np.meshgrid(
        np.linspace(35, 55, 6), np.linspace(-10, 20, 8), indexing="ij"
    )
    ext_locs = np.column_stack(
        [lat_g.ravel(), lon_g.ravel(), np.full(lat_g.size, R_E + 110e3)]
    )
    int_locs = np.column_stack(
        [lat_g.ravel(), lon_g.ravel(), np.full(lat_g.size, R_E - 500e3)]
    )

    # Synthetic truth: internal currents mirror a fraction of the external
    ext_amps = rng.normal(0, 1e4, ext_locs.shape[0])
    int_amps = -0.35 * ext_amps

    nobs = 30
    obs_locs = np.column_stack(
        [rng.uniform(36, 54, nobs), rng.uniform(-9, 19, nobs), np.full(nobs, R_E)]
    )
    B_obs = np.tensordot(ext_amps, T_df(obs_locs, ext_locs), (0, 2)) + np.tensordot(
        int_amps, T_df(obs_locs, int_locs), (0, 2)
    )
    B_obs += rng.normal(0, 0.02 * np.max(np.abs(B_obs)), B_obs.shape)

    # Fit both shells at once by stacking their SEC locations
    both = SECS(sec_df_loc=np.vstack([ext_locs, int_locs]))
    both.fit(obs_locs, B_obs, epsilon=0.01)

    # Split the fitted amplitudes back onto individual shells for prediction
    n_ext = len(ext_locs)
    ext_fit = SECS(sec_df_loc=ext_locs)
    ext_fit.sec_amps = both.sec_amps[:, :n_ext]
    int_fit = SECS(sec_df_loc=int_locs)
    int_fit.sec_amps = both.sec_amps[:, n_ext:]

    # Compare along a latitude profile
    profile = np.column_stack(
        [np.linspace(38, 52, 100), np.full(100, 5.0), np.full(100, R_E)]
    )
    B_ext_true = np.tensordot(ext_amps, T_df(profile, ext_locs), (0, 2))
    B_int_true = np.tensordot(int_amps, T_df(profile, int_locs), (0, 2))

    B_ext_pred = ext_fit.predict_B(profile)
    B_int_pred = int_fit.predict_B(profile)

    fig, ax = plt.subplots()
    lats = profile[:, 0]
    ax.plot(lats, B_ext_true[:, 0] * 1e9, "C0", label="external truth")
    ax.plot(lats, B_ext_pred[:, 0] * 1e9, "C0--", label="external fit")
    ax.plot(lats, B_int_true[:, 0] * 1e9, "C1", label="internal truth")
    ax.plot(lats, B_int_pred[:, 0] * 1e9, "C1--", label="internal fit")
    ax.set_xlabel("latitude (deg)")
    ax.set_ylabel(r"B$_x$ (nT)")
    ax.set_title("Two-shell separation of external and internal sources")
    ax.legend()
    plt.show()


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.121 seconds)


.. _sphx_glr_download_examples_plot_internal_external.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_internal_external.ipynb <plot_internal_external.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_internal_external.py <plot_internal_external.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_internal_external.zip <plot_internal_external.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_