{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Getting Started with PsiData\n\nThis example is intended as a quick reference for interacting with the :func:`~psi_io.mhd_io.PsiData`\nlazy reader. The nuances of this API's design and full capabilities is explored in the subsequent\nexamples.\n\n.. seealso::\n   `sphx_glr_gallery_03_mhd_io_p02_psidata_reading.py` for the full\n   read/vslice API, and\n   `sphx_glr_gallery_01_reading_files_p01_reading_file_meta_data.py`\n   for metadata inspection.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import warnings\n\nfrom psi_io import PsiData\nfrom psi_data import fetch_mas_data, fetch_pot3d_data\n\nfrom psi_io.mhd_io import MetaDataWarning"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Fetching Data\n\nFor the purposes of this example, the :func:`~psi_data.fetch_mas_data` and\n:func:`~psi_data.fetch_pot3d_data` functions, found in the :mod:`psi_data` package, are\nused to retrieve example MAS and POT3D HDF files. These functions return a simple namespace of\nfilepaths for each variable.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mas_files= fetch_mas_data(domains='cor', variables='br,vr')\npot3d_file = fetch_pot3d_data(variables='br')\n\nprint(f\"Example MAS Br file        \u2192    .../{mas_files.cor_br.name}\")\nprint(f\"Example MAS Vr file        \u2192    .../{mas_files.cor_vr.name}\")\nprint(f\"Example POT3D Br file      \u2192    .../{pot3d_file.br.name}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Creating a Reader\n\nThe :func:`~psi_io.mhd_io.PsiData` reader API requires a filepath and ``model`` name; if the\n``model`` name is not recognized, *i.e.* not ``mas`` or ``pot3d``, metadata cannot be inferred\nand must be supplied explicitly (either as keyword arguments to :func:`~psi_io.mhd_io.PsiData`\nor as attributes written to the HDF dataset itself - see :func:`~psi_io.psi_io.write_hdf_meta`\nfor more information on how to write metadata to HDF datasets.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mas_br = PsiData(mas_files.cor_br, model='mas')\nprint(f\"{mas_br!r}\")\n\nwith warnings.catch_warnings(record=True) as w:\n    pot3d_br = PsiData(pot3d_file.br, model='pot3d')\n    print(f\"METADATA WARNING: {w[-1].message}\")\n\nprint(f\"{pot3d_br!r}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-danger\"><h4>Warning</h4><p>POT3D applies no normalization to its output. The stored values are in whatever physical\n   unit the input photospheric magnetogram used \u2014 most commonly Gauss, but this is not encoded\n   in the file. The fallback unit for POT3D is dimensionless_unscaled (scale factor 1).\n   As a result, it is incumbent upon the user to pass in the correct unit at instantiation\n   time, or to set it manually after the fact, to ensure accurate physical interpretation and\n   unit-aware computations.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pot3d_br.unit = 'Gauss'  # manually set the unit to Gauss for accurate physical interpretation\nprint(f\"Updated POT3D Br unit: {pot3d_br.unit}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Passing the Wrong Model or Mesh\n\n:func:`~psi_io.mhd_io.PsiData` does **not** verify that the declared ``model`` actually\nmatches the file \u2014 it trusts ``model`` and infers every piece of metadata (quantity name and\ndescription, unit normalization, and mesh staggering) from that model's mapping together with\nthe filename. Declaring the wrong model therefore *silently* attaches the wrong metadata.\n\nBelow we deliberately open the POT3D ``br`` file as a MAS quantity. The reader happily\nconstructs, but the inferred unit is now ``MAS_b`` instead of ``POT3D_b`` \u2014 a different\nnormalization scale factor that would corrupt any physical-unit computation \u2014 and the mesh\nstagger is taken from MAS' ``br`` descriptor rather than POT3D's:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "with warnings.catch_warnings(record=True) as w:\n    warnings.simplefilter(\"always\")\n    err_reader = PsiData(pot3d_file.br, model='mas')\n\nprint(f\"Inferred unit (wrong) : {err_reader.unit}               (should be POT3D_b)\")\nprint(f\"Inferred mesh (wrong) : {err_reader.mesh}    (should be MAIN, HALF, HALF)\")\nfor warning in w:\n    if issubclass(warning.category, MetaDataWarning):\n        print(f\"METADATA WARNING: {warning.message}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-danger\"><h4>Warning</h4><p>**A model/mesh mismatch is mostly silent.** The wrong *unit* is applied with no warning at\n   all, so values look numerically fine but are scaled incorrectly. The only automatic signal\n   is a best-effort :class:`~psi_io.mhd_io.MetaDataWarning` raised during *scale* validation: a\n   ``t`` or ``p`` coordinate declared on the *main* mesh but holding non-zero values at the\n   inner boundary carries the tell-tale offset of a *half*-mesh axis, which is what trips the\n   warnings above. This heuristic catches some \u2014 but not all \u2014 staggering mistakes, and it\n   never inspects the data unit.\n\n   The mesh code is not cosmetic: it drives the cell-averaging performed by\n   ``read(mesh=...)`` and the half-/main-mesh offsets used to bracket coordinates in\n   :meth:`~psi_io.mhd_io.PsiData.vslice`. A wrong code averages along the wrong axis or shifts\n   the coordinate\u2013data alignment by up to half a grid cell, again without raising an error.\n\n   When the model cannot be inferred \u2014 or is wrong \u2014 supply the correct metadata explicitly via\n   the ``unit`` and ``mesh`` keywords at instantiation (or set them afterward, as with the unit\n   above), or write the correct attributes into the HDF dataset with\n   :func:`~psi_io.psi_io.write_hdf_meta`.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Putting It Together: Reading on the Main Mesh in Physical Units\n\nThe most common end-to-end task is simply *read the whole field, on the main mesh, in physical\nunits*. Passing ``mesh='main'`` averages any half-mesh axis onto the main mesh, and\n``unit='physical'`` decomposes the data to CGS base units. :meth:`~psi_io.mhd_io.PsiData.read`\nreturns the data array together with the (remeshed) ``r``, ``t``, ``p`` coordinate scales.\n\nWe also use the reader as a **context manager**. ``PsiData`` holds an open OS-level handle to\nthe underlying HDF file (opened when the reader is constructed). Using it in a ``with`` block\nguarantees that handle is closed *deterministically* the moment the block exits \u2014 including if\nan exception is raised mid-read \u2014 rather than lingering until the object is garbage-collected.\nThis matters when looping over many files, where leaked handles can otherwise exhaust the\nprocess's file-descriptor limit or keep files locked. ``__enter__`` (re)opens the handle and\nreturns the reader; ``__exit__`` closes it.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "with PsiData(mas_files.cor_vr, model='mas') as reader:\n    data, r, t, p = reader.read(mesh='main', unit='physical')\n    print(f\"quantity   : {reader.name} [{reader.desc}]\")\n    print(f\"data       : shape={data.shape}, unit={data.unit}\")\n    print(f\"r / t / p  : {r.shape[0]} \u00d7 {t.shape[0]} \u00d7 {p.shape[0]} points\")\n\n# Outside the ``with`` block the underlying HDF handle has been closed, but the data we read is\n# an ordinary in-memory array and remains fully usable:\nprint(f\"max |Vr|   : {abs(data).max():.3e}\")"
      ]
    }
  ],
  "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.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}