{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# MHDweb Integration Part II\n\nThis is the second of a two-part series.\n`sphx_glr_gallery_99_advanced_plots_p09_integrating_mhdweb_p1.py`\nqueried the [MHDweb REST API](https://predsci.com/mhdweb2/api) to download\ncoronal and heliospheric $(B_r, B_\\theta, B_\\phi)$ field files and the\nSolar Orbiter spacecraft connectivity mapping; this example loads those\noutputs and produces four visualizations:\n\n1. A radially scaled overview of $B_r$ in both model domains.\n2. A static scene of the coronal $B_r$ slice with spacecraft positions,\n   ballistically mapped positions, and backward-traced coronal connectivity.\n3. An animation of inter-domain and spacecraft mapping traces from a fixed\n   wide-angle heliospheric perspective.\n4. A close-up coronal animation whose camera follows each spacecraft's\n   longitude through the sequence.\n\nMagnetic connectivity is computed in two ways:\n\n- **Spacecraft mapping** \u2014 :class:`~mapflpy.tracer.TracerMP` traces backward\n  from the ballistically mapped positions at\n  $r_1 \\approx 30\\,R_\\odot$ through the coronal domain.  The\n  spacecraft position is then prepended to form a continuous path from\n  the heliosphere to the solar surface.\n- **Inter-domain tracing** \u2014 :func:`~mapflpy.scripts._inter_domain_tracing`\n  launches field-line integration from the spacecraft positions directly,\n  crossing the coronal\u2013heliospheric domain boundary to produce end-to-end\n  connectivity from the spacecraft to $r_0 = 1\\,R_\\odot$.\n\n.. seealso::\n\n   `sphx_glr_gallery_99_advanced_plots_p09_integrating_mhdweb_p1.py`\n      Part I \u2014 queries MHDweb, downloads field files, and saves the\n      spacecraft mapping table.\n\n   `sphx_glr_gallery_02_stack_mesh_mixin_p03_fieldlines.py`\n      Introduction to fieldline rendering with\n      :meth:`~pyvisual.core.mixins.StackMeshMixin.add_fieldlines`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os\nfrom contextlib import ExitStack\nfrom math import pi\nfrom pathlib import Path\nfrom zipfile import ZipFile\nfrom mapflpy.tracer import TracerMP\nfrom mapflpy.scripts import _inter_domain_tracing\nfrom mapflpy.utils import combine_and_pad_fieldlines\nimport astropy.units as u\nfrom sunpy.sun.constants import sidereal_rotation_rate\nfrom astropy.table import QTable\nimport numpy as np\n\nfrom pyvisual import Plot3d, SphericalMesh"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Extract and Load Field Data\n\nThe ZIP archives downloaded in Part I are extracted to their respective\ndirectories. :class:`~pyvisual.core.mesh3d.SphericalMesh` loads each\n``br002.h5`` file directly.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "OUTPUT_DIR = Path(os.environ.get(\"STATIC_ASSETS\", \"\")).resolve()\nCOR_OUTPUT_DIR = OUTPUT_DIR / 'cor_mag_field'\nHEL_OUTPUT_DIR = OUTPUT_DIR / 'hel_mag_field'\n\nprint(\"Extracting coronal magnetic field files...\")\nwith ZipFile(COR_OUTPUT_DIR / 'cor_mag_field.zip') as cor_zip:\n    print(f'cor_mag_field.zip namelist: {cor_zip.namelist()}')\n    cor_zip.extractall(path=COR_OUTPUT_DIR)\nprint(\"Extracting heliospheric magnetic field files...\")\nwith ZipFile(HEL_OUTPUT_DIR / 'hel_mag_field.zip') as hel_zip:\n    print(f'hel_mag_field.zip namelist: {hel_zip.namelist()}')\n    hel_zip.extractall(path=OUTPUT_DIR / 'hel_mag_field')\n\ncor_br_mesh = SphericalMesh(COR_OUTPUT_DIR / 'br002.h5')\nhel_br_mesh = SphericalMesh(HEL_OUTPUT_DIR / 'br002.h5')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Parse the Spacecraft Mapping Table\n\nThe ECSV file written in Part I is read into an\n:class:`astropy.table.QTable`, which preserves column units (distances in\n$R_\\odot$, angles in radians).\n\nThe three position arrays \u2014 each of shape $(3, N)$ in\n$(r, \\theta, \\phi)$ order \u2014 correspond to the three legs of the\nspacecraft connectivity chain returned by the\n[MHDweb spacecraft-mapping endpoint](https://predsci.com/mhdweb2/api):\n\n- ``spacecraft_positions`` \u2014 actual spacecraft location.\n- ``balmapped_positions`` \u2014 position ballistically mapped to the outer\n  coronal boundary $r_1 \\approx 30\\,R_\\odot$.\n- ``traced_positions`` \u2014 magnetic footpoint at the inner boundary\n  $r_0 = 1\\,R_\\odot$.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "spacecraft_mapping = QTable.read(OUTPUT_DIR / 'spacecraft_mapping.ecsv', format='ascii.ecsv')\n\nspacecraft_positions = np.stack(tuple(spacecraft_mapping[f'sc_pos_{dim}'].value for dim in 'rtp'))\nbalmapped_positions = np.stack(tuple(spacecraft_mapping[f'r1_pos_{dim}'].value for dim in 'rtp'))\ntraced_positions = np.stack(tuple(spacecraft_mapping[f'r0_pos_{dim}'].value for dim in 'rtp'))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To properly visualize the ballistic mappings, we need to construct a continuous\npath from the spacecraft position through the heliosphere to the coronal footpoint.\n\nFor each time step, a radial path is constructed from 50 linearly spaced radial\npositions (beginning from the spacecraft position and ending at $30\\,R_\\odot$).\nThe time shift for each radial position is computed based on the in situ flow speed, and the\ncorresponding longitudinal shift is calculated using\n:mod:`sunpy`'s :class:`~sunpy.sun.constants.sidereal_rotation_rate`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "balmapping_radial_path = np.linspace(spacecraft_mapping['sc_pos_r'].value, 30, 50) * u.R_sun\ntime_shift = (spacecraft_mapping['sc_pos_r'] - balmapping_radial_path) / (spacecraft_mapping['flow_speed'])\nlongitudinal_shift = time_shift * sidereal_rotation_rate\nbalmapping_longitudinal_path = ((spacecraft_mapping['sc_pos_p'] + longitudinal_shift) % (360 * u.deg)).to(u.rad)\n\nballistic_mapping_trajectory = np.stack(\n    (balmapping_radial_path.value,\n     np.full_like(balmapping_radial_path.value, spacecraft_positions[1]),\n     balmapping_longitudinal_path.value),\n    axis=1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Trace Magnetic Connectivity\n\n:class:`~mapflpy.tracer.TracerMP` is a multiprocessing-capable tracer that\nmust be used as a context manager; :class:`contextlib.ExitStack` manages\ntwo tracer contexts simultaneously so that both are cleanly shut down even\nif an exception occurs.\n\n**Spacecraft mapping traces** \u2014 backward integration from\n``balmapped_positions`` (at $r_1$) through the coronal domain.  The\nspacecraft positions are prepended along axis 0 so the resulting array\nrepresents the complete path from the heliosphere through to the coronal\nfootpoint.\n\n**Inter-domain traces** \u2014 :func:`~mapflpy.scripts._inter_domain_tracing`\nlaunches from the actual spacecraft positions and integrates through both\nthe heliospheric and coronal domains, crossing the domain boundary\nautomatically.  :func:`~mapflpy.utils.combine_and_pad_fieldlines` merges\nthe per-domain segments into a single NaN-padded array suitable for\n:meth:`~pyvisual.core.mixins.StackMeshMixin.add_splines`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "with ExitStack() as cstack:\n    cor_tracer = cstack.enter_context(\n        TracerMP(*[COR_OUTPUT_DIR / f'b{dim}002.h5' for dim in 'rtp'], context='fork'))\n    hel_tracer = cstack.enter_context(\n        TracerMP(*[HEL_OUTPUT_DIR / f'b{dim}002.h5' for dim in 'rtp'], context='fork'))\n\n    spacecraft_mapping_traces = cor_tracer.trace_bwd(\n        launch_points=balmapped_positions)\n    spacecraft_mapping_traces = np.concatenate(\n        (ballistic_mapping_trajectory, spacecraft_mapping_traces.geometry))\n\n    inter_domain_traces, *_ = _inter_domain_tracing(cor_tracer, hel_tracer,\n        launch_points=spacecraft_positions)\n    inter_domain_traces = combine_and_pad_fieldlines(inter_domain_traces)\n\ncommon_kwargs = dict(cmap='seismic', clim=10, show_scalar_bar=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Two-Domain $B_r$ Overview\n\nBoth coronal and heliospheric $B_r$ meshes are rendered in a single scene\nto provide context for the connectivity visualizations that follow.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plotter = Plot3d()\nplotter.add_axes()\nplotter.add_mesh(cor_br_mesh.radially_scale(), opacity=0.8, **common_kwargs)\nplotter.add_mesh(hel_br_mesh.radially_scale(), opacity=0.6, **common_kwargs)\nplotter.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Spacecraft Connectivity in the Corona\n\nSpacecraft positions (colored by time index) and\nballistically mapped positions are shown as point clouds; the backward\ntraces connecting them through the coronal domain are rendered as splines.\n:func:`numpy.moveaxis` transposes the geometry array from\n$(M, 3, N)$ to $(3, M, N)$ for unpacking into\n:meth:`~pyvisual.core.mixins.StackMeshMixin.add_splines`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plotter = Plot3d()\nplotter.add_axes()\nplotter.add_mesh(cor_br_mesh.slice(normal='x', origin=(1, 0, 0)), opacity=0.8, **common_kwargs)\nplotter.add_points(*spacecraft_positions,\n                   np.arange(len(spacecraft_mapping)),\n                   point_size=3)\nplotter.add_points(*balmapped_positions,\n                   np.arange(len(spacecraft_mapping)),\n                   point_size=3)\nplotter.add_splines(*np.moveaxis(spacecraft_mapping_traces, 1, 0),\n                    np.arange(len(spacecraft_mapping)))\nplotter.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Animate Spacecraft Mapping in the Heliosphere\n\nA fixed wide-field observer at\n$(r, \\theta, \\phi) = (400\\,R_\\odot,\\;\\pi/4,\\;0)$ with a\n$200^\\circ$ field of view frames the entire heliospheric domain.\nEach frame advances one time step, drawing the spacecraft mapping trace\n(white) and inter-domain trace (red) in a rolling buffer of five named\nactors so the scene never accumulates more than five trace pairs.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plotter = Plot3d()\nplotter.add_axes()\nplotter.add_points(*spacecraft_positions,\n                   np.arange(len(spacecraft_mapping)),\n                   point_size=3)\nplotter.add_mesh(cor_br_mesh[1, ...], **common_kwargs)\nplotter.observer_position = 400, pi / 4, 0\nplotter.observer_fov_view = 200\nplotter.open_gif(OUTPUT_DIR / 'spacecraft_mapping_heliosphere.gif', fps=40)\nfor i, (scmap_trace, interdomain_trace) in enumerate(\n        zip(spacecraft_mapping_traces.T, inter_domain_traces.T)):\n    plotter.add_spline(*scmap_trace, name=f'scmap_trace_{i % 5}', color='white', line_width=3)\n    plotter.add_spline(*interdomain_trace, name=f'interdomain_trace_{i % 5}', color='red', line_width=3)\n    plotter.write_frame()\nplotter.close()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Animate Spacecraft Mapping in the Corona\n\nThe same traces are re-animated from a close-up coronal perspective.\nThe observer is placed at $r = 15\\,R_\\odot$ near the ecliptic\n($\\theta = 3\\pi/8$) and its longitude tracks each spacecraft\nposition plus a $\\pi/6$ offset, keeping the active trace near the\ncenter of the $4^\\circ$ field of view as the sequence progresses.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plotter = Plot3d()\nplotter.add_axes()\nplotter.add_points(*spacecraft_positions,\n                   np.arange(len(spacecraft_mapping)),\n                   point_size=3)\nplotter.add_mesh(cor_br_mesh[1, ...], **common_kwargs)\nplotter.open_gif(OUTPUT_DIR / 'spacecraft_mapping_corona.gif', fps=40)\nfor i, (scmap_trace, interdomain_trace) in enumerate(\n        zip(spacecraft_mapping_traces.T, inter_domain_traces.T)):\n    plotter.observer_position = 15, 3 * pi / 8, spacecraft_positions[2, i] + pi / 6\n    plotter.observer_fov_view = 4\n    plotter.add_spline(*scmap_trace, name=f'scmap_trace_{i % 5}', color='white', line_width=3)\n    plotter.add_spline(*interdomain_trace, name=f'interdomain_trace_{i % 5}', color='red', line_width=3)\n    plotter.write_frame()\nplotter.close()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.13.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}