{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Polarity-Inversion-Line Seeded Fieldlines\n\nThis example traces magnetic fieldlines seeded from the polarity inversion\nline (PIL) \u2014 the curve on the solar surface where the radial magnetic field\n$B_r$ changes sign.  The PIL marks the boundary between open positive\nand negative flux and is a natural seed surface for displaying the global\ntopology of the coronal field. It can be used to locate the apex of coronal\narcades and/or highlight bald-patch topologies.\n\nThe workflow exploits a key property of\n:class:`~pyvisual.core.mesh3d.SphericalMesh`: because it wraps a\n:class:`pyvista.RectilinearGrid` whose internal axes store\n$(r, \\theta, \\phi)$ directly, PyVista's\n:meth:`pyvista.DataSetFilters.contour` returns a\n:class:`pyvista.PolyData` whose ``points`` array is already in\n$(r, \\theta, \\phi)$ coordinates.  This means the contour points can\nbe passed straight to\n:func:`~mapflpy.scripts.run_fwdbwd_tracing` as ``launch_points`` without\nany intermediate coordinate conversion.\n\nSee also `sphx_glr_gallery_02_stack_mesh_mixin_p03_fieldlines.py` for a\ngeneral introduction to fieldline rendering, and\n`sphx_glr_gallery_99_advanced_plots_p01_combining_multiple_elements.py` for a\nbroader scene that layers slices, contours, and fieldlines.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom mapflpy.scripts import run_fwdbwd_tracing\nfrom pyvisual import Plot3d\nfrom pyvisual.core.mesh3d import SphericalMesh\nfrom pyvisual.utils.data import fetch_datasets"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Load Data and Extract the Polarity Inversion Line\n\nBuild a full :class:`~pyvisual.core.mesh3d.SphericalMesh` from the coronal\n$B_r$ HDF file.  Index ``mesh[5, ...]`` selects a 2-D spherical shell\nat radial index 5 \u2014 a thin :class:`pyvista.RectilinearGrid` that spans the\nfull $(\\\\theta, \\\\phi)$ domain at a fixed $r$.  The ellipsis\n(``...``) expands to fill the remaining two axes.\n\n:meth:`pyvista.DataSetFilters.contour` with ``isosurfaces=[0]`` then finds\nthe iso-curve where $B_r = 0$ on that shell.  The result is the\npolarity inversion line as a :class:`pyvista.PolyData` of line segments\nwhose ``points`` array holds $(r, \\\\theta, \\\\phi)$ coordinates \u2014\nbecause the internal axes of :class:`~pyvisual.core.mesh3d.SphericalMesh`\n*are* the spherical coordinates, not Cartesian positions.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mag_field = fetch_datasets(\"cor\", [\"br\", \"bt\", \"bp\"])\n\nmesh = SphericalMesh(mag_field.cor_br)\nneutraline = mesh[5, ...].contour(isosurfaces=[0])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Trace Fieldlines from the PIL\n\n:func:`~mapflpy.scripts.run_fwdbwd_tracing` integrates the magnetic field\nin both the forward and backward directions from each seed point, so that\nevery fieldline has footpoints on both the inner\n($r = 1\\,R_{\\\\odot}$) and outer ($r = 30\\,R_{\\\\odot}$) boundaries.\nThis bidirectional tracing is required for accurate polarity classification\nand ensures that closed fieldlines connecting two inner-boundary footpoints\nare also captured.\n\n``launch_points=neutraline.points.T`` passes the PIL vertices directly as\nthe $(3, N)$ seed array in $(r, \\\\theta, \\\\phi)$ order.  No\ncoordinate conversion is needed because the contour inherits the\n:class:`~pyvisual.core.mesh3d.SphericalMesh` coordinate convention.\n\n:func:`numpy.moveaxis` transposes the ``(M, 3, N)`` geometry array to\n$(3, M, N)$ so that unpacking with ``*`` feeds the three coordinate\ncomponents directly to :meth:`~pyvisual.core.mixins.StackMeshMixin.add_fieldlines`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "traces = run_fwdbwd_tracing(*mag_field, launch_points=neutraline.points.T, context='fork')\nr, t, p = np.moveaxis(traces.geometry, 1, 0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Render the Scene\n\nThree layers are combined in a single :class:`~pyvisual.core.plot3d.Plot3d`\nscene:\n\n- The radial shell ``mesh[5, ...]`` coloured by $B_r$, providing\n  context for the polarity structure at the seed radius.\n- The PIL rendered as a white tube, showing the exact seed curve.\n- The fieldline bundle, each line assigned a random hue via\n  ``coloring='random'``, illustrating the global connectivity of the coronal\n  field anchored at the polarity boundary.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plotter = Plot3d()\nplotter.show_axes()\nplotter.add_sun()\nplotter.add_mesh(mesh[5, ...], cmap='seismic', clim=(-1e-1, 1e-1), opacity=0.5, show_scalar_bar=False)\nplotter.add_mesh(neutraline, color='white', line_width=3, render_lines_as_tubes=True)\nplotter.add_fieldlines(r, t, p, coloring='random', line_width=1, show_scalar_bar=False)\nplotter.observer_focus = 0, 0, 0\nplotter.observer_fov_view = 10\nplotter.show()"
      ]
    }
  ],
  "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
}