diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md new file mode 100644 index 0000000..0703ccf --- /dev/null +++ b/.github/CONTRIBUTING.md @@ -0,0 +1,93 @@ +# Contributing Guidelines + +Thank you for your interest in the project. Everything is still a work in progress and +your assistance/feedback will inform future development! + +## Opening an issue + +If you find a bug, want a new feature, or anything else, please create an issue. +Remember to: + +- **Be descriptive.** + If you are reporting a bug, include instructions on how to replicate it. + If you had an installation error, include the exact steps you were trying & details on your system. +- **Search for relevant issues.** + Look through the open (and recently closed) issues to see if others have reported something similar. + +## Testing + +`glowpython2` is tested on Linux, MacOS and Windows. +Currently three Python versions are used (stable, latest, and experimental). + +The tests are run as a GitHub action, where the source file & instructions are located in [.github/workflows/ci.yml](workflows/ci.yml) + +### Installation for running tests + +Because `glowpython2` uses meson-python as its build backend, +**editable installs (`pip install -e .`) do not work correctly for running the test suite**. +Data files required by the Fortran backend are not placed where the package expects them. +Use a regular install instead: + +```bash +pip install .[test] +``` + +**After modifying Fortran source**, you must reinstall for changes to take effect: + +```bash +pip install . +``` + +### Running tests locally + +From the repo root: + +```bash +pytest tests/ +``` + +You may run an individual test file if the tests are taking too long. + +### Test layout + +- **`tests/test_examples.py`** — automatically runs every script in `Examples/` as a subprocess smoke test. No maintenance needed; new examples are picked up automatically. +- **`tests/test_no_precip.py`**, **`test_maxwellian.py`** — structure and physics assertions for the top-level GLOW functions. +- **`tests/test_iri.py`** — assertions for both `Iri90` and `Iri2020`. +- **`tests/test_msis.py`** — assertions for both `NrlMsis00` and `NrlMsis21`. +- **`tests/conftest.py`** — shared fixtures (`sample_time`, `sample_lat`, `sample_lon`). + +### Adding or modifying tests + +- **New example script**: drop a `.py` file in `Examples/` — it will be picked up by `test_examples.py` automatically. +- **New assertion test**: add a `test_*.py` file in `tests/`, or add a test function to an existing file. Use the fixtures from `conftest.py` for the standard time/location inputs. +- **Physics assertions** should check that outputs are physically reasonable (positive densities, sane temperature ranges, expected species present) rather than matching exact numerical values, so tests stay valid across minor model updates. + +## Changelog + +All notable changes should be recorded in `CHANGELOG.md` before a release is tagged. + +**During development:** add entries under the `[Unreleased]` section at the top of `CHANGELOG.md` as you work. Use the appropriate subsection: + +- **Added** — new features or capabilities +- **Changed** — changes to existing behavior +- **Fixed** — bug fixes +- **Removed** — removed features or APIs +- **Deprecated** — features that will be removed in a future release + +**At release time:** +1. Rename `[Unreleased]` to `[x.y.z] - YYYY-MM-DD` matching the version tag you're about to create. +2. Add a new empty `[Unreleased]` section above it. +3. Tag the release in git: `git tag vx.y.z`. + +**Example entry:** + +```markdown +## [Unreleased] +### Added +- Support for custom altitude grids + +### Fixed +- Incorrect electron density below 80 km under certain solar conditions +``` + +Do not list every commit — only changes a user of the library would care about. Trivial internal refactors, test updates, and CI fixes can be omitted. diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..8f6c38c --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,71 @@ +name: CI + +on: + push: + branches: [master, develop] + pull_request: + branches: [master] + +jobs: + test: + name: Run tests (${{ matrix.os }}, Python ${{ matrix.python-version }}) + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, macos-latest, windows-latest] + python-version: ["3.11", "3.12", "3.14"] + + steps: + - uses: actions/checkout@v4 + + - name: Set up Python (Linux/macOS) + if: runner.os != 'Windows' + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + + - name: Install gfortran (Linux) + if: runner.os == 'Linux' + run: sudo apt-get install -y gfortran + + - name: Install gfortran (macOS) + if: runner.os == 'macOS' + run: brew install gcc + + - name: Set up conda with Python (Windows) + if: runner.os == 'Windows' + uses: conda-incubator/setup-miniconda@v3 + with: + miniforge-version: latest + python-version: ${{ matrix.python-version }} + + - name: Create short temp directory (Windows) + if: runner.os == 'Windows' + run: New-Item -ItemType Directory -Force -Path C:\T | Out-Null + + - name: Install gfortran + glowpython2 (Windows) + if: runner.os == 'Windows' + env: + TEMP: C:\T + TMP: C:\T + run: | + conda install -c conda-forge gfortran -y + conda run --no-capture-output pip install .[test] + conda run python -c "import glowpython2; print(glowpython2.__version__)" + + - name: Install (Linux/macOS) + if: runner.os != 'Windows' + run: pip install .[test] + + - name: Tests (Linux/macOS) + if: runner.os != 'Windows' + run: python -m pytest tests + + - name: Tests (Windows) + # Running all tests on Windows for now... + # If the runtime starts getting too long, this can be changed to just one test + if: runner.os == 'Windows' + run: conda run --no-capture-output python -m pytest tests + # run: conda run --no-capture-output python -m pytest tests/test_no_precip.py + # ^^ If we just want to run one test diff --git a/.gitignore b/.gitignore index b34ac17..a383f28 100644 --- a/.gitignore +++ b/.gitignore @@ -35,3 +35,4 @@ OLDGLOW/ # Others *.txt *.png +!Examples/*.png diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..9276c8b --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,33 @@ +# Changelog + +All notable changes to glowpython2 are documented here. + +## [Unreleased] +### Added +- `CHANGELOG.md` with history migrated from README +- GitHub Actions CI workflow for Linux, macOS, and Windows +- Functionality test suite with `pytest` + +### Changed +- `CONTRIBUTING.md` expanded with changelog workflow and guidelines +- Changelog section removed from `README.md` +- README install instructions updated for Windows and general clarity +- `python-dateutil` dependency removed; minimum Python version bumped +- `pyproject.toml` dependencies tidied +- Moved files from tests/ to Examples & updates README to use local versions of some images + +### Fixed +- Windows install: + - gfortran now recommended to be installed via conda instead of MSYS2 + - pip forced to use shorter paths to avoid path-length issues + +## [0.0.3] - 2026-01-09 +### Changed +- Dataset structure changed (**breaking**) + +### Fixed +- Version bumped to 0.0.3 + +## [0.0.1] - 2025-11-01 +### Added +- Initial release diff --git a/Examples/IriComparison.py b/Examples/IriComparison.py new file mode 100644 index 0000000..48da809 --- /dev/null +++ b/Examples/IriComparison.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python +""" +Compare IRI-1990 and IRI-2020 ion densities and temperatures. +""" +import matplotlib +import matplotlib.pyplot as plt +from datetime import datetime, UTC + +from glowpython2 import Iri90 +from glowpython2.utils import alt_grid +from iri20py import Iri2020 + +matplotlib.rcParams.update({'mathtext.fontset': 'cm'}) +matplotlib.rc('font', **{'family': 'serif', 'serif': ['Times New Roman']}) + +time = datetime(2015, 12, 13, 10, 0, 0, tzinfo=UTC) +lat = 42.6 +lon = -71.2 +alt = alt_grid() + +_, ds90 = Iri90().evaluate(time=time, lat=lat, lon=lon, alt=alt) +_, ds20 = Iri2020().evaluate(time=time, lat=lat, lon=lon, alt=alt) + +fig, ax = plt.subplots(figsize=(8, 6), dpi=300) +tax = ax.twiny() + +species = ['O+', 'O2+', 'N+', 'NO+'] +descs = ['O^+', 'O_2^+', 'N^+', 'NO^+'] +colors = ['r', 'g', 'b', 'm'] +lines, labels = [], [] + +for spec, color, desc in zip(species, colors, descs): + l20, = ds20[spec].plot(y='alt_km', ax=ax, color=color) + l90, = ds90[spec].plot(y='alt_km', ax=ax, linestyle='--', color=color) + lines += [l20, l90] + labels += [f'IRI-20 ${desc}$', f'IRI-90 ${desc}$'] + +for ds, ls, label_te, label_ti in [ + (ds20, '-', 'IRI-20 $T_e$', 'IRI-20 $T_i$'), + (ds90, '--', 'IRI-90 $T_e$', 'IRI-90 $T_i$'), +]: + l_te, = ds['Te'].plot(y='alt_km', ax=tax, color='k', linestyle=ls) + l_ti, = ds['Ti'].plot(y='alt_km', ax=tax, color='c', alpha=0.7, linestyle=ls) + lines += [l_te, l_ti] + labels += [label_te, label_ti] + +ax.set_title('IRI-90 vs IRI-20') +ax.set_xscale('log') +ax.set_xlabel('Number Density [cm$^{-3}$]') +tax.set_xlabel('Temperature [K]') +ax.set_xlim(1e-3, None) +tax.set_xlim(100, None) +ax.legend(lines, labels, loc='upper left', fontsize='small') + +plt.savefig('iri_comparison.png', dpi=300, bbox_inches='tight') +plt.show() diff --git a/Examples/ModelComparison.py b/Examples/ModelComparison.py new file mode 100644 index 0000000..d0a4691 --- /dev/null +++ b/Examples/ModelComparison.py @@ -0,0 +1,141 @@ +#!/usr/bin/env python +""" +Compare GLOW model output across atmospheric model configurations: + - MSIS-2000 + IRI-1990 + - MSIS-2.1 + IRI-2020 + - MSIS-2.1 + IRI-2020 + ModGLOW coefficients +""" +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from datetime import datetime +from matplotlib.gridspec import GridSpec + +from glowpython2 import no_precipitation + +matplotlib.rcParams.update({'mathtext.fontset': 'cm'}) +matplotlib.rc('font', **{'family': 'serif', 'serif': ['Times New Roman']}) + +time = datetime(2022, 3, 22, 22, 0, 0) +glat = 42.6 +glon = -71.2 +Nbins = 250 + +versions = ['MSIS00_IRI90', 'MSIS21_IRI20', 'MODGLOW'] +biglabels = ['MSIS-2000 + IRI-1990', 'MSIS-2.1 + IRI-2020', 'MSIS-2.1 + IRI-2020 + ModGLOW'] +lstyles = ['-', '--', '-.'] +alphas = [0.5, 0.7, 0.5] + +ionos = {} +for version in versions: + if version == 'MODGLOW': + iono = no_precipitation(time, glat, glon, Nbins, version='MSIS21_IRI20', magmodel='IGRF14', newcoeffs=True) + else: + magmodel = 'IGRF14' if version == 'MSIS21_IRI20' else 'POGO68' + iono = no_precipitation(time, glat, glon, Nbins, version=version, magmodel=magmodel) + ionos[version] = iono + +# Density plot +grid = GridSpec(2, 2, height_ratios=[0.05, 1], width_ratios=[1, 1], hspace=0.1, wspace=0.15) +den_figure = plt.figure(figsize=(6.4, 4.8), dpi=300) +den_axs = np.asarray([den_figure.add_subplot(grid[1, 0]), den_figure.add_subplot(grid[1, 1])]) +legend_ax = den_figure.add_subplot(grid[0, :]) +legend_ax.axis('off') + +descs = { + "O": "O", "O2": "O$_2$", "N2": "N$_2$", "NO": "NO", + "O+": "O$^+$", "O2+": "O$_2^+$", "NO+": "NO$^+$", "N(2D)": "N($^2$D)", 'NeIn': 'e$^-$' +} +den_lines = [{}, {}] +for version, ls, alpha in zip(versions, lstyles, alphas): + iono = ionos[version] + den_axs[0].set_prop_cycle(None) + for v in ("O", "O2", "N2", "NO"): + l, = den_axs[0].plot(iono[v], iono[v].alt_km, linestyle=ls, lw=0.75, alpha=alpha) + if v not in den_lines[0]: + den_lines[0][v] = l + den_axs[1].set_prop_cycle(None) + for v in ("O+", "O2+", "NO+", "N(2D)", 'NeIn'): + l, = den_axs[1].plot(iono[v], iono[v].alt_km, linestyle=ls, lw=0.75, alpha=alpha) + if v not in den_lines[1]: + den_lines[1][v] = l + +for ax, lines, title in zip(den_axs, den_lines, ['Neutrals', 'Ions']): + ax.set_xscale('log') + ax.set_ylim(60, 1000) + ax.set_xlim(1, None) + ax.grid(True) + ax.set_title(title, fontsize='medium') + ax.legend([v for v in lines.values()], [descs.get(k) for k in lines], loc='best', fontsize='small') +den_axs[0].set_ylabel('Altitude [km]') +den_axs[1].yaxis.set_ticklabels([]) +for lab, ls, alpha in zip(biglabels, lstyles, alphas): + legend_ax.plot([], [], label=lab, color='k', linestyle=ls, alpha=alpha) +legend_ax.legend(loc='center', ncol=3, fontsize='small', mode='expand', frameon=False) +den_figure.text(0.5, 0.04, 'Number Density [cm$^{-3}$]', ha='center') +den_figure.savefig('model_comparison_density.png', dpi=300, bbox_inches='tight') + +# Temperature plot +temp_grid = GridSpec(2, 1, height_ratios=[0.05, 1], hspace=0.1) +temp_figure = plt.figure(figsize=(6.4, 4.8), dpi=300) +temp_ax = temp_figure.add_subplot(temp_grid[1, 0]) +temp_legend_ax = temp_figure.add_subplot(temp_grid[0, 0]) +temp_legend_ax.axis('off') +temp_lines = {} +for version, ls, alpha in zip(versions, lstyles, alphas): + iono = ionos[version] + temp_ax.set_prop_cycle(None) + for v in ('Te', 'Ti', 'Tn'): + l, = temp_ax.plot(iono[v], iono[v].alt_km, linestyle=ls, lw=0.75, alpha=alpha) + if v not in temp_lines: + temp_lines[v] = l +temp_ax.set_xlabel('Temperature [K]') +temp_ax.set_ylabel('Altitude [km]') +temp_ax.set_title('Temperature') +temp_ax.grid(True) +temp_ax.set_ylim(60, 1000) +temp_ax.set_xlim(100, None) +temp_ax.legend(temp_lines.values(), [f'${k}$' for k in temp_lines], loc='best') +for lab, ls, alpha in zip(biglabels, lstyles, alphas): + temp_legend_ax.plot([], [], label=lab, color='k', linestyle=ls, alpha=alpha) +temp_legend_ax.legend(loc='center', ncol=3, fontsize='small', mode='expand', frameon=False) +temp_figure.savefig('model_comparison_temperature.png', dpi=300, bbox_inches='tight') + +# VER plot +ver_grid = GridSpec(2, 3, height_ratios=[0.05, 1], hspace=0.1, wspace=0.15) +ver_figure = plt.figure(figsize=(6.4, 4.8), dpi=300) +ver_axs = [ver_figure.add_subplot(ver_grid[1, i]) for i in range(3)] +ver_legend_ax = ver_figure.add_subplot(ver_grid[0, :]) +ver_legend_ax.axis('off') +ver_components = { + 'Visible': ["4278", "5577", "6300", "5200"], + 'Infrared': ["7320", "7774", "8446", "10400"], + 'Ultraviolet': ["1304", "1356", "1493", "3371", "3644", "3726", "LBH"], +} +ver_lines = {k: {} for k in ver_components} +for version, ls, alpha in zip(versions, lstyles, alphas): + iono = ionos[version] + for ax, kind in zip(ver_axs, ver_components): + ax.set_prop_cycle(None) + for line in ver_components[kind]: + l, = ax.plot(iono['ver'].sel(wavelength=line), iono['ver'].alt_km, linestyle=ls, lw=0.75, alpha=alpha) + if line not in ver_lines[kind]: + ver_lines[kind][line] = l +for ax, kind in zip(ver_axs, ver_components): + ax.set_title(kind, fontsize='medium') + ax.grid(True) + ax.set_ylim(60, 1000) + ax.set_xlim(1e-5, None) + ax.set_xscale('log') + labels = [f'{k} Å' if k != 'LBH' else 'LBH' for k in ver_lines[kind]] + ax.legend(ver_lines[kind].values(), labels, loc='best', fontsize='small') +ver_axs[0].set_ylabel('Altitude [km]') +for ax in ver_axs[1:]: + ax.yaxis.set_ticklabels([]) +for lab, ls, alpha in zip(biglabels, lstyles, alphas): + ver_legend_ax.plot([], [], label=lab, color='k', linestyle=ls, alpha=alpha) +ver_legend_ax.legend(loc='center', ncol=3, fontsize='small', mode='expand', frameon=False) +ver_figure.text(0.5, 0.04, 'Volume Emission Rate [cm$^{-3}$ s$^{-1}$]', ha='center') +ver_figure.savefig('model_comparison_ver.png', dpi=300, bbox_inches='tight') + +plt.show() diff --git a/Examples/MsisComparison.py b/Examples/MsisComparison.py new file mode 100644 index 0000000..e0f4d96 --- /dev/null +++ b/Examples/MsisComparison.py @@ -0,0 +1,50 @@ +#!/usr/bin/env python +""" +Compare NRLMSISE-00 and NRLMSIS-2.1 neutral densities and temperatures. +""" +import matplotlib +import matplotlib.pyplot as plt +from datetime import datetime, UTC + +from glowpython2 import NrlMsis00 +from glowpython2.utils import alt_grid +from msis21py import NrlMsis21 + +matplotlib.rcParams.update({'mathtext.fontset': 'cm'}) +matplotlib.rc('font', **{'family': 'serif', 'serif': ['Times New Roman']}) + +time = datetime(2022, 1, 13, 10, 0, 0, tzinfo=UTC) +lat = 0.0 +lon = 0.0 +alt = alt_grid() + +ds00 = NrlMsis00().evaluate(time=time, lat=lat, lon=lon, alt=alt) +ds21 = NrlMsis21().evaluate(time=time, lat=lat, lon=lon, alt=alt) + +fig, ax = plt.subplots(figsize=(8, 6), dpi=300) +tax = ax.twiny() + +species = ['O', 'O2', 'N2', 'NO'] +descs = ['O', 'O$_2$', 'N$_2$', 'NO'] +colors = ['r', 'g', 'b', 'm'] +lines, labels = [], [] + +for spec, color, desc in zip(species, colors, descs): + l21, = ds21[spec].plot(y='alt_km', ax=ax, color=color) + l00, = ds00[spec].plot(y='alt_km', ax=ax, linestyle='--', color=color) + lines += [l21, l00] + labels += [f'MSIS-21 {desc}', f'MSIS-00 {desc}'] + +l21, = ds21['Tn'].plot(y='alt_km', ax=tax, color='k') +l00, = ds00['Tn'].plot(y='alt_km', ax=tax, color='k', linestyle='--') +lines += [l21, l00] +labels += ['MSIS-21 $T_n$', 'MSIS-00 $T_n$'] + +ax.set_title('NRLMSISE-00 vs NRLMSIS-2.1') +ax.set_xscale('log') +ax.set_xlabel('Number Density [cm$^{-3}$]') +tax.set_xlabel('Neutral Temperature [K]') +ax.legend(lines, labels) + +plt.savefig('msis_comparison.png', dpi=300, bbox_inches='tight') +plt.show() diff --git a/Examples/NoPrecipitation.py b/Examples/NoPrecipitation.py index 735f9d7..d8af49c 100644 --- a/Examples/NoPrecipitation.py +++ b/Examples/NoPrecipitation.py @@ -17,8 +17,8 @@ iono = glow.no_precipitation(time, glat, glon, Nbins) # %% simple plots if False: - plot.precip(iono["precip"]) # all zeros as intended - plot.altitude(iono) + Plot.precip(iono["precip"]) # all zeros as intended + Plot.altitude(iono) plot = Plot() plot.density(iono) diff --git a/Examples/iri_comparison.png b/Examples/iri_comparison.png new file mode 100644 index 0000000..334198f Binary files /dev/null and b/Examples/iri_comparison.png differ diff --git a/Examples/model_comparison_density.png b/Examples/model_comparison_density.png new file mode 100644 index 0000000..e395866 Binary files /dev/null and b/Examples/model_comparison_density.png differ diff --git a/Examples/model_comparison_temperature.png b/Examples/model_comparison_temperature.png new file mode 100644 index 0000000..f993359 Binary files /dev/null and b/Examples/model_comparison_temperature.png differ diff --git a/Examples/model_comparison_ver.png b/Examples/model_comparison_ver.png new file mode 100644 index 0000000..0718046 Binary files /dev/null and b/Examples/model_comparison_ver.png differ diff --git a/Examples/msis_comparison.png b/Examples/msis_comparison.png new file mode 100644 index 0000000..f4f9251 Binary files /dev/null and b/Examples/msis_comparison.png differ diff --git a/README.md b/README.md index 83ac2eb..f394d8e 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,11 @@ [![DOI](https://zenodo.org/badge/1026267765.svg)](https://zenodo.org/badge/latestdoi/1026267765) # GlowPython2 -The FORTRAN GLobal airglOW ([GLOW](https://github.com/NCAR/GLOW)) Model in Python ≥ 3.9. +The FORTRAN GLobal airglOW ([GLOW](https://github.com/NCAR/GLOW)) Model in Python. A Fortran compiler is **REQUIRED**. -Note: This version uses `meson` and `ninja` as the build system, and does not rely on `distutils`, -and is Python 3.12 compatible. +Note: This version uses `meson` and `ninja` as the build system, does not rely on `distutils`, +and is compatible with Python ≥ 3.11. This library also allows parallelization of model evaluation using the [multiprocessing](https://docs.python.org/3/library/multiprocessing.html) module. @@ -22,55 +22,116 @@ into the following fundamental steps: ## Prerequisites A Fortran compiler is **REQUIRED**. + ### Linux Ensure that you have the following development packages installed: - `build-essential` (for `gcc`, `g++`, `make`, etc.) - `gfortran` (Fortran compiler) + ### macOS Ensure that you have the Xcode Command Line Tools installed. You can install them by running: ```sh xcode-select --install ``` + Install [`homebrew`](https://brew.sh/) if you haven't already, and then install `gfortran`: ```sh brew install gfortran ``` + Note: For macOS Big Sur and above, you may need to add the following line to your environment script (`~/.zshrc` if using ZSH, or the relevant shell init script): + ```sh export LIBRARY_PATH="$LIBRARY_PATH:/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/lib" ``` Then reopen the terminal. This fixes the issue where `-lSystem` fails for `gfortran`. - + ### Windows (amd64 or x86_64 targets) -On Windows, MSYS2 is the preferred distribution for installing the Fortran compiler toolchain (for GNU Compiler Collection). -- Install [MSYS2](https://www.msys2.org/). - Note the directory where MSYS2 was installed (defaults to `C:\msys64`) [referred to as `MSYS_INSTALL_DIR`]. - It is not recommended to change this directory. -- Launch the MSYS2 terminal (MSYS2 MSYS application on the start menu) -- Update MSYS2 environment (assuming fresh install): - ```sh - pacman -Syu # Restart the terminal - pacman -Su # Update packages - ``` -- Install the GNU Compiler Collection: - ```sh - pacman -S --needed base-devel mingw-w64-ucrt-x86_64-toolchain mingw-w64-ucrt-x86_64-gcc-fortran - ``` -- Add `MSYS_INSTALL_DIR\ucrt64\bin` (defaults to `C:\msys64\ucrt64\bin`) to `PATH`: - - Search for `env` in the Start menu, - - Select "Edit the system environment variables", - - Click "Environment Variables", - - Double click 'Path' under 'User variables for USER' - - Click "New" - - Type in, or paste the full path to `ucrt64\bin` (defaults to `C:\msys64\ucrt64\bin`) - - Click "Ok" on the environment variable windows to save the changes. -- Continue with installation instructions for the Python packages below, in a new, regular terminal (e.g. Command Prompt or PowerShell with Python installed). + +> [!IMPORTANT] +> **Windows Path Length Limitation**: `glowpython2` and its dependencies (`msis21py`, +> `iri20py`) compile Fortran code using Meson & pip, which can put temporary object +> files with very long paths. Windows has a 260-character path limit by default. +> +> You **must set the TEMP and TMP environment variables to short paths** before +> installing, regardless of which method you use to install `gfortran`. + +**Conda (Recommended):** + +Using [conda](https://docs.conda.io/en/latest/) or [miniforge](https://github.com/conda-forge/miniforge): + +```sh +conda install -c conda-forge gfortran +mkdir C:\T +set TEMP=C:\T +set TMP=C:\T +``` + +**Other Options:** + +
+MSYS2 + +1. Install [MSYS2](https://www.msys2.org/) (defaults to `C:\msys64`) +2. Launch the MSYS2 terminal from the start menu +3. Update MSYS2: + ```sh + pacman -Syu # Restart the terminal when prompted + pacman -Su # Update packages + ``` +4. Install the GNU Compiler Collection: + ```sh + pacman -S --needed base-devel mingw-w64-ucrt-x86_64-toolchain mingw-w64-ucrt-x86_64-gcc-fortran + ``` +5. Add MSYS2 to PATH: + - Search for "env" in the Start menu + - Select "Edit the system environment variables" + - Click "Environment Variables" + - Double-click 'Path' under 'User variables' + - Click "New" + - Add `C:\msys64\ucrt64\bin` + - Click OK to save +6. In a new regular terminal (Command Prompt or PowerShell), set short temp paths and install: + ```sh + mkdir C:\T + set TEMP=C:\T + set TMP=C:\T + ``` + +
+ +
+Winlibs + +1. Download a pre-built MinGW-w64 bundle from [winlibs.com](https://winlibs.com) +2. Extract it anywhere on your system +3. Add the `bin\` folder to your PATH: + - Search for "env" in the Start menu + - Select "Edit the system environment variables" + - Click "Environment Variables" + - Double-click 'Path' under 'User variables' + - Click "New" + - Add the full path to the `bin` folder (e.g., `C:\mingw\bin`) + - Click OK to save +4. In a new terminal, set short temp paths and install: + ```sh + mkdir C:\T + set TEMP=C:\T + set TMP=C:\T + ``` + +
+ + +Then continue with the instructions below. Either `pip install` or use build manually. > [!NOTE] -> Change the toolchain names accordingly for Windows arm64. -> This platform has not been tested and is not officially supported. +> Windows arm64 is not tested and is not officially supported. ## Installation + +**Please ensure you have followed the prequisite steps above before installing!** + Direct installation using pip: ```sh pip install glowpython2 @@ -78,22 +139,31 @@ pip install glowpython2 Install from source repository: ```sh -pip install glowpython2@git+https://github.com/sunipkm/glowpython2 +git clone https://github.com/sunipkm/glowpython2 +cd glowpython2 +pip install . ``` +> **Note:** editable installs (`pip install -e .`) are not supported — meson-python's +> editable mode does not correctly stage the data files required by the Fortran backend. +> After modifying Fortran source, re-run `pip install .` to rebuild. +> +> If you just want to download the latest commit & build it from source, use: + +> ```sh +> pip install glowpython2@git+https://github.com/sunipkm/glowpython2 +> ``` + Requires (and installs) [geomagdata](https://pypi.org/project/geomagdata/) for timezone aware geomagnetic parameter retrieval. + ## Usage + ### Command-line examples + * `Glow2Maxwellian`: asdf * `Glow2NoPrecip`: asdf -### Pre-defined examples - -* [`Maxwellian.py`](https://raw.githubusercontent.com/sunipkm/glowpython2/refs/heads/master/Examples/Maxwellian.py): Maxwellian precipitation, specify Q (flux) and E0 (characteristic energy). -* [`NoPrecipitation.py`](https://raw.githubusercontent.com/sunipkm/glowpython2/refs/heads/master/Examples/NoPrecipitation.py): No precipitating electrons. -* [`test_glow.py`](https://raw.githubusercontent.com/sunipkm/glowpython2/refs/heads/master/tests/test_glow.py): Compares GLOW results between `glowpython` and `glowpython2` modules without precipitation, as well as with NRLMSISE-2.1 atmosphere and IRI-2020 ionosphere. - -These examples are also available on the command line: `Glow2Maxwellian` and `Glow2NoPrecip`. +> This is a work in progress! ### Using `glowpython2` module In Python source code, import the module and call the `maxwellian` function: @@ -105,20 +175,29 @@ iono = glow.maxwellian(time, glat, glon, Nbins, Q, Echar) Read the module documentation for more information. -### Example Plots +### Examples + + +Example scripts are in [Examples/](Examples/). +These examples cover no-precipitation and Maxwellian runs, GLOW model comparisons across atmosphere/ionosphere configurations, and side-by-side IRI and MSIS model comparisons. + | Densities | Temperatures | | :---: | :---: | -| ![Comparison of [`glowpython`](https://pypi.org/project/glowpython/) and `glowpython2` densities with no precipitation.](https://raw.githubusercontent.com/sunipkm/glowpython2/refs/heads/master/tests/glow_den.png) | ![Comparison of [`glowpython`](https://pypi.org/project/glowpython/) and `glowpython2` temperatures with no precipitation.](https://raw.githubusercontent.com/sunipkm/glowpython2/refs/heads/master/tests/glow_temp.png) | +| ![Model comparison: neutral densities.](Examples/model_comparison_density.png) | ![Model comparison: temperatures.](Examples/model_comparison_temperature.png) | + +| Volume Emission Rates | IRI Comparison | MSIS Comparison | +| :---: | :---: | :---: | +| ![Model comparison: volume emission rates.](Examples/model_comparison_ver.png) | ![IRI-90 vs IRI-2020 ion densities and temperatures.](Examples/iri_comparison.png) | ![MSIS-00 vs MSIS-2.1 neutral densities.](Examples/msis_comparison.png) | -| Volume Emission Rates | -| :---: | -| ![Comparison of [`glowpython`](https://pypi.org/project/glowpython/) and `glowpython2` volume emission rates with no precipitation.](https://raw.githubusercontent.com/sunipkm/glowpython2/refs/heads/master/tests/glow_ver.png) | ### Options + GlowPython2 has a few configuration options at runtime - - Magnetic field model: Use `POGO68` for GLOW 0.98 behavior, `IGRF14` for the [IGRF-14](https://www.ncei.noaa.gov/products/international-geomagnetic-reference-field) model (extracted from [IRI-2020](https://irimodel.org/)). - Atmosphere and ionosphere model: Use `MSIS00_IRI90` for GLOW 0.98 behavior, `MSIS21_IRI20` for [NRLMSIS 2.1](https://map.nrl.navy.mil/map/pub/nrl/NRLMSIS/NRLMSIS2.1/) and [IRI-2020](https://irimodel.org/). - New coefficients from [ModGLOW implementation](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2017JA025026): Set `newcoeffs` to `True`. This updates the reaction rates with data from newer publications. + ### Model Output + The returned output is a [xarray.Dataset](http://xarray.pydata.org/en/stable/generated/xarray.Dataset.html) containing outputs from GLOW: @@ -181,11 +260,13 @@ containing outputs from GLOW: All available keys carry unit and description. # Acknowledgement + Sincerest gratitude to Dr. Stan Solomon for developing, and open-sourcing the GLOW model. We are very thankful to Guy Gribbs for providing the [`gchem_modglow.f90`](https://raw.githubusercontent.com/sunipkm/glowpython2/refs/heads/master/src/GLOW/gchem_modglow.f90) file with the modified reaction rate coefficients for the `GCHEM` subroutine. His prompt response saved us many hours of work. # Citation + If you use this code in your work, please cite the repository: ```bibtex @software{sunipkm_glowpython2_2025, @@ -199,8 +280,3 @@ If you use this code in your work, please cite the repository: url = {https://github.com/sunipkm/glowpython2}, } ``` - -# Changelog -- 2026-01-09 - - Version bumped to 0.0.3 - - BREAKING: Dataset structure changed. diff --git a/pyproject.toml b/pyproject.toml index 3514e2b..e849279 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,7 +13,7 @@ build-backend = "mesonpy" [project] name = "glowpython2" dynamic = ["version"] -requires-python = ">=3.9" +requires-python = ">=3.11" authors = [ {name = "Sunip K. Mukherjee", email = "sunipkmukherjee@gmail.com"}, ] @@ -34,16 +34,18 @@ classifiers = [ license = {file = "LICENSE"} readme = "README.md" dependencies = [ - 'meson', - 'ninja', - 'numpy >= 1.14.5', - 'xarray >= 0.15', - 'geomagdata >= 1.6.1', + 'numpy >= 1.14.5', + 'xarray >= 0.15', + 'geomagdata >= 1.6.1', 'pytz', - 'msis21py>=0.0.5', - 'iri20py>=0.0.5', + 'matplotlib', + 'msis21py >= 0.0.5', + 'iri20py >= 0.0.5', ] +[project.optional-dependencies] +test = ["pytest"] + [tool.setuptools_scm] [project.urls] diff --git a/src/glowpython2/examples.py b/src/glowpython2/examples.py index 71fd2ac..b891395 100644 --- a/src/glowpython2/examples.py +++ b/src/glowpython2/examples.py @@ -4,7 +4,7 @@ from .base import maxwellian, no_precipitation from .plots import Plot from matplotlib.pyplot import show -from dateutil.parser import parse +from datetime import datetime import sys import argparse import json @@ -27,7 +27,7 @@ def Maxwellian(): args = parser.parse_args() - time = parse(args.time) + time = datetime.fromisoformat(args.time) glat = args.glat glon = args.glon # %% flux [erg cm-2 s-1 == mW m-2 s-1] @@ -75,7 +75,7 @@ def NoPrecipitation(): args = parser.parse_args() - time = parse(args.time) + time = datetime.fromisoformat(args.time) glat = args.glat glon = args.glon Nbins = args.Nbins diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..63fa450 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,23 @@ +import os +import pytest +from datetime import datetime, UTC + +os.environ.setdefault("MPLBACKEND", "Agg") + + +@pytest.fixture(scope="session") +def sample_time(): + # Define a single time for all tests + return datetime(2015, 12, 13, 10, 0, 0, tzinfo=UTC) + + +@pytest.fixture(scope="session") +def sample_lat(): + # Single latitude for tests + return 65.1 + + +@pytest.fixture(scope="session") +def sample_lon(): + # Longitude for all tests + return -147.5 diff --git a/tests/glow_den.png b/tests/glow_den.png deleted file mode 100644 index 7472731..0000000 Binary files a/tests/glow_den.png and /dev/null differ diff --git a/tests/glow_temp.png b/tests/glow_temp.png deleted file mode 100644 index ac12aa3..0000000 Binary files a/tests/glow_temp.png and /dev/null differ diff --git a/tests/glow_ver.png b/tests/glow_ver.png deleted file mode 100644 index 3557958..0000000 Binary files a/tests/glow_ver.png and /dev/null differ diff --git a/tests/test_examples.py b/tests/test_examples.py new file mode 100644 index 0000000..032f0dd --- /dev/null +++ b/tests/test_examples.py @@ -0,0 +1,23 @@ +import os +import subprocess +import sys +from pathlib import Path + +import pytest + +EXAMPLES_DIR = Path(__file__).parent.parent / "Examples" +EXAMPLES = sorted(EXAMPLES_DIR.glob("*.py")) + +# Runs all of the python files in Examples/ +# Just a smoke test, does not actually verify their outputs. +# Just makes sure they can be run! +# New files detected auromatically +@pytest.mark.parametrize("script", EXAMPLES, ids=[e.name for e in EXAMPLES]) +def test_example(script): + result = subprocess.run( + [sys.executable, str(script)], + env=os.environ, + capture_output=True, + cwd=script.parent, + ) + assert result.returncode == 0, result.stderr.decode() diff --git a/tests/test_glow.py b/tests/test_glow.py deleted file mode 100644 index 7c12e31..0000000 --- a/tests/test_glow.py +++ /dev/null @@ -1,225 +0,0 @@ -# %% -from datetime import UTC -from matplotlib.gridspec import GridSpec -import pytz -from glowpython import no_precipitation as glownoprecip -from glowpython2 import no_precipitation -import matplotlib -import matplotlib.pyplot as plt -import numpy as np -from dateutil.parser import parse -# %% -usetex = False -if not usetex: - # computer modern math text - matplotlib.rcParams.update({'mathtext.fontset': 'cm'}) - -matplotlib.rc( - 'font', **{ - 'family': 'serif', - 'serif': ['Times' if usetex else 'Times New Roman'] - } -) -# for Palatino and other serif fonts use: -# rc('font',**{'family':'serif','serif':['Palatino']}) -matplotlib.rc('text', usetex=usetex) -# %% -time = parse('2022-3-22T22:00:00') -glat = 42.6 -glon = -71.2 -Nbins = 250 -tec = None -versions = ['MSIS00_IRI90', 'MSIS21_IRI20', 'MODGLOW'] -biglabels = ['MSIS-2000 + IRI-1990', 'MSIS-2.1 + IRI-2020', 'MSIS-2.1 + IRI-2020 + ModGLOW'] -lstyles = ['-', '--', '-.'] -alphas = [0.5, 0.7, 0.5] -ionos = {} -for version in versions: - if version == 'GLOW': - iono = glownoprecip(time, glat, glon, Nbins) - elif version == 'MODGLOW': - iono = no_precipitation( - time, glat, glon, Nbins, tec=tec, - version='MSIS21_IRI20', magmodel='IGRF14', - newcoeffs=True - ) - else: - magmodel = 'IGRF14' if version == 'MSIS21_IRI20' else 'POGO68' - iono = no_precipitation(time, glat, glon, Nbins, tec=tec, version=version, magmodel=magmodel) # type: ignore - ionos[version] = iono -# %% Density plot -grid = GridSpec(2, 2, height_ratios=[0.05, 1], width_ratios=[1, 1], hspace=0.1, wspace=0.15) -den_figure = plt.figure(figsize=(6.4, 4.8), dpi=300) -den_axs = np.asarray([den_figure.add_subplot(grid[1, 0]), den_figure.add_subplot(grid[1, 1])]) -legend_axs = den_figure.add_subplot(grid[0, :]) -legend_axs.axis('off') -den_lines = [{}, {}] # neutrals, ions -descs = { - "O": "O", "O2": "O$_2$", "N2": "N$_2$", "NO": "NO", - "O+": "O$^+$", "O2+": "O$_2^+$", "NO+": "NO$^+$", "N(2D)": "N($^2$D)", 'NeIn': 'e$^-$' -} -for version, linestyle, alpha in zip(versions, lstyles, alphas): - iono = ionos[version] - ax = den_axs[0] - ax.set_prop_cycle(None) # reset color cycle - for v in ("O", "O2", "N2", "NO"): - l, = ax.plot(iono[v], iono[v].alt_km, label=f'{version} {v}', linestyle=linestyle, lw=0.75, alpha=alpha) - if v not in den_lines[0]: - den_lines[0][v] = l - -for version, linestyle, alpha in zip(versions, lstyles, alphas): - iono = ionos[version] - ax = den_axs[1] - ax.set_prop_cycle(None) # reset color cycle - for v in ("O+", "O2+", "NO+", "N(2D)", 'NeIn'): - l, = ax.plot(iono[v], iono[v].alt_km, label=f'{version} {v}', linestyle=linestyle, lw=0.75, alpha=alpha) - if v not in den_lines[1]: - den_lines[1][v] = l - -den_axs[0].set_xscale("log") -den_axs[0].set_ylim(60, 1000) -den_axs[0].set_ylabel("Altitude [km]") -den_axs[0].set_title("Neutrals", fontsize='medium') -den_axs[0].grid(True) -den_axs[0].set_xlim(1, None) -lines = [] -labels = [] -for k, v in den_lines[0].items(): - lines.append(v) - labels.append(descs.get(k)) -den_axs[0].legend(lines, labels, loc="best", fontsize='small') -lines = [] -labels = [] -den_axs[1].set_xscale("log") -den_axs[1].set_title("Ions", fontsize='medium') -den_axs[1].grid(True) -den_axs[1].set_xlim(1, None) -den_axs[1].yaxis.set_ticklabels([]) -for k, v in den_lines[1].items(): - lines.append(v) - labels.append(descs.get(k)) -den_axs[1].legend(lines, labels, loc="best", fontsize='small') -for lab, ls, alpha in zip(biglabels, lstyles, alphas): - legend_axs.plot([], [], label=lab, color='k', linestyle=ls, alpha=alpha) -legend_axs.legend(loc='center', ncol=3, fontsize='small', mode='expand', frameon=False) -den_figure.text(0.5, 0.04, "Number Density [cm$^{-3}$]", ha='center') -den_figure.text(0.5, 0.9, time.astimezone(UTC).astimezone(pytz.timezone('US/Eastern')).isoformat(sep=' '), ha='center') -den_figure.savefig('glow_den.png', dpi=300, bbox_inches='tight') -plt.show() -# %% -# Temperature plot -temp_grid = GridSpec(2, 1, height_ratios=[0.05, 1], hspace=0.1) -temp_figure = plt.figure(figsize=(6.4, 4.8), dpi=300) -temp_ax = temp_figure.add_subplot(temp_grid[1, 0]) -temp_legend_ax = temp_figure.add_subplot(temp_grid[0, 0]) -temp_legend_ax.axis('off') -temp_lines = {} -for version, linestyle, alpha in zip(versions, lstyles, alphas): - iono = ionos[version] - ax = temp_ax - ax.set_prop_cycle(None) # reset color cycle - l, = ax.plot(iono['Te'], iono['Te'].alt_km, label=f'{version} Te', linestyle=linestyle, lw=0.75, alpha=alpha) - if 'Te' not in temp_lines: - temp_lines['Te'] = l - l, = ax.plot(iono['Ti'], iono['Ti'].alt_km, label=f'{version} Ti', linestyle=linestyle, lw=0.75, alpha=alpha) - if 'Ti' not in temp_lines: - temp_lines['Ti'] = l - l, = ax.plot(iono['Tn'], iono['Tn'].alt_km, label=f'{version} Tn', linestyle=linestyle, lw=0.75, alpha=alpha) - if 'Tn' not in temp_lines: - temp_lines['Tn'] = l -temp_ax.set_xlabel("Temperature [K]") -temp_ax.set_ylabel("altitude [km]") -temp_ax.set_title("Temperature") -temp_ax.grid(True) -temp_ax.set_ylim(60, 1000) -temp_ax.set_xlim(100, None) -lines = [] -labels = [] -for k, v in temp_lines.items(): - lines.append(v) - labels.append(f'${k}$') -temp_ax.legend(lines, labels, loc="best") -for lab, ls, alpha in zip(biglabels, lstyles, alphas): - temp_legend_ax.plot([], [], label=lab, color='k', linestyle=ls, alpha=alpha) -temp_legend_ax.legend(loc='center', ncol=3, fontsize='small', mode='expand', frameon=False) -temp_figure.text(0.5, 0.9, time.astimezone(UTC).astimezone(pytz.timezone('US/Eastern')).isoformat(sep=' '), ha='center') -temp_figure.savefig('glow_temp.png', dpi=300, bbox_inches='tight') -plt.show() -# %% VER plot -ver_grid = GridSpec(2, 3, height_ratios=[0.05, 1], hspace=0.1, wspace=0.15) -ver_figure = plt.figure(figsize=(6.4, 4.8), dpi=300) -ver_axs = np.asarray([ - ver_figure.add_subplot(ver_grid[1, 0]), - ver_figure.add_subplot(ver_grid[1, 1]), - ver_figure.add_subplot(ver_grid[1, 2]) -]) -ver_legend_ax = ver_figure.add_subplot(ver_grid[0, :]) -ver_legend_ax.axis('off') -ver_lines = { - 'Visible': {}, - 'Infrared': {}, - 'Ultraviolet': {} -} -ver_components = { - 'Visible': ["4278", "5577", "6300", "5200"], - 'Infrared': ["7320", "7774", "8446", "10400"], - 'Ultraviolet': ["1304", "1356", "1493", "3371", "3644", "3726", "LBH"] -} -for version, linestyle, alpha in zip(versions, lstyles, alphas): - iono = ionos[version] - for ax, kind in zip(ver_axs, ('Visible', 'Infrared', 'Ultraviolet')): - ax.set_prop_cycle(None) # reset color cycle - for line in ver_components[kind]: - l, = ax.plot(iono['ver'].sel(wavelength=line), iono['ver'].alt_km, - label=f'{version} {line} Å', linestyle=linestyle, lw=0.75, alpha=alpha) - if line not in ver_lines[kind]: - ver_lines[kind][line] = l - -for ax, kind in zip(ver_axs, ('Visible', 'Infrared', 'Ultraviolet')): - ax.set_title(f"{kind}", fontsize='medium') - ax.grid(True) - ax.set_ylim(60, 1000) - ax.set_xlim(1e-5, None) - ax.set_xscale("log") - lines = [] - labels = [] - for k, v in ver_lines[kind].items(): - lines.append(v) - if k != 'LBH': - labels.append(f'{k} Å') - else: - labels.append('LBH') - ax.legend(lines, labels, loc="best", fontsize='small') -for ax in ver_axs[1:]: - ax.yaxis.set_ticklabels([]) -ver_axs[0].set_ylabel("Altitude [km]") -ver_figure.text(0.5, 0.04, "Volume Emission Rate [cm$^{-3}$ s$^{-1}$]", ha='center') -ver_figure.text(0.5, 0.9, time.astimezone(UTC).astimezone(pytz.timezone('US/Eastern')).isoformat(sep=' '), ha='center') -for lab, ls, alpha in zip(biglabels, lstyles, alphas): - ver_legend_ax.plot([], [], label=lab, color='k', linestyle=ls, alpha=alpha) -ver_legend_ax.legend(loc='center', ncol=3, fontsize='small', mode='expand', frameon=False) -ver_figure.savefig('glow_ver.png', dpi=300, bbox_inches='tight') -plt.show() - -# %% -iono = ionos['MSIS00_IRI90'] -ver_above = iono['ver'].sel(alt_km=slice(300, None)) -ver_below = iono['ver'].sel(alt_km=slice(None, 300)) -alt_above = ver_above.alt_km.values -alt_below = ver_below.alt_km.values -# nan filter -loc = np.where(np.isnan(ver_above.values)) -ver_above.values[loc] = 0.0 -loc = np.where(np.isnan(ver_below.values)) -ver_below.values[loc] = 0.0 -# ver_above_integrated = np.trapezoid(ver_above.values, x=(alt_above*1e2), axis=0) -# ver_below_integrated = np.trapezoid(ver_below.values, x=(alt_below*1e2), axis=0) -ver_above = ver_above.integrate(coord='alt_km') -# ver_above.values = ver_above_integrated -ver_below = ver_below.integrate(coord='alt_km') -# ver_below.values = ver_below_integrated -for wavelength in ['5577', '6300']: - print(f'Wavelength: {wavelength} Å') - print(f' VER above 300 km: {ver_above.sel(wavelength=wavelength).values:.3e} cm^-2 s^-1') - print(f' VER below 300 km: {ver_below.sel(wavelength=wavelength).values:.3e} cm^-2 s^-1') -# %% diff --git a/tests/test_iri.py b/tests/test_iri.py new file mode 100644 index 0000000..0d02a9d --- /dev/null +++ b/tests/test_iri.py @@ -0,0 +1,115 @@ +import numpy as np +import pytest +from glowpython2 import Iri90 +from glowpython2.utils import alt_grid +from iri20py import Iri2020 + +SPECIES = ["Ne", "O+", "O2+", "NO+"] +TEMPERATURES = ["Te", "Ti"] + + +@pytest.fixture(scope="module") +def ds90(sample_time, sample_lat, sample_lon): + _, ds = Iri90().evaluate(time=sample_time, lat=sample_lat, lon=sample_lon, alt=alt_grid()) + return ds + + +@pytest.fixture(scope="module") +def ds20(sample_time, sample_lat, sample_lon): + _, ds = Iri2020().evaluate(time=sample_time, lat=sample_lat, lon=sample_lon, alt=alt_grid()) + return ds + + +# --- output structure --- + +@pytest.mark.parametrize("ds_fixture", ["ds90", "ds20"]) +def test_has_alt_coord(ds_fixture, request): + ds = request.getfixturevalue(ds_fixture) + assert "alt_km" in ds.coords + + +@pytest.mark.parametrize("ds_fixture,species", [ + (f, s) for f in ["ds90", "ds20"] for s in SPECIES +]) +def test_has_species(ds_fixture, species, request): + ds = request.getfixturevalue(ds_fixture) + assert species in ds, f"{ds_fixture} missing {species}" + + +@pytest.mark.parametrize("ds_fixture,temp", [ + (f, t) for f in ["ds90", "ds20"] for t in TEMPERATURES +]) +def test_has_temperature(ds_fixture, temp, request): + ds = request.getfixturevalue(ds_fixture) + assert temp in ds, f"{ds_fixture} missing {temp}" + + +# --- physical sanity --- + +@pytest.mark.parametrize("ds_fixture", ["ds90", "ds20"]) +def test_electron_density_valid_range(ds_fixture, request): + # Ne must be positive in the F-region + ds = request.getfixturevalue(ds_fixture) + ne = ds["Ne"].sel(alt_km=slice(150, 500)).values + assert np.all(ne > 0) + + +@pytest.mark.parametrize("ds_fixture", ["ds90", "ds20"]) +def test_electron_density_fill_below_boundary(ds_fixture, request): + # Below ~80 km IRI is undefined; must be exactly the -1e-6 fill + ds = request.getfixturevalue(ds_fixture) + ne_low = ds["Ne"].sel(alt_km=slice(None, 79)).values + assert np.all(ne_low == pytest.approx(-1e-6, rel=1e-3)) + + +@pytest.mark.parametrize("ds_fixture", ["ds90", "ds20"]) +def test_temperatures_valid_range(ds_fixture, request): + # Te and Ti must be physically reasonable above 150 km + ds = request.getfixturevalue(ds_fixture) + for temp in TEMPERATURES: + vals = ds[temp].sel(alt_km=slice(150, 500)).values + assert np.all(vals > 100), f"{ds_fixture} {temp} unreasonably cold above 150 km" + + +def test_temperatures_fill_below_boundary_iri90(ds90): + # IRI-90 uses -1.0 as fill below 120 km; must be exactly that + for temp in TEMPERATURES: + vals = ds90[temp].sel(alt_km=slice(None, 119)).values + assert np.all(vals == -1.0), f"ds90 {temp} unexpected values below IRI lower boundary" + + +@pytest.mark.parametrize("ds_fixture", ["ds90", "ds20"]) +def test_alt_range(ds_fixture, request): + ds = request.getfixturevalue(ds_fixture) + alt = ds.coords["alt_km"].values + assert alt.min() >= 50 + assert alt.max() <= 2000 + + +@pytest.mark.parametrize("ds_fixture", ["ds90", "ds20"]) +def test_te_ti_physically_reasonable(ds_fixture, request): + # Te and Ti above 200 km should be warm but not absurd + ds = request.getfixturevalue(ds_fixture) + for temp in TEMPERATURES: + vals = ds[temp].sel(alt_km=slice(200, 500)).values + assert np.all(vals > 500), f"{ds_fixture} {temp} unreasonably cold above 200 km" + assert np.all(vals < 15000), f"{ds_fixture} {temp} unreasonably hot above 200 km" + + +@pytest.mark.parametrize("ds_fixture", ["ds90", "ds20"]) +def test_ne_peaks_in_f_region(ds_fixture, request): + # NmF2 must fall in the F-region; anything outside 150-600 km is a model failure + ds = request.getfixturevalue(ds_fixture) + ne = ds["Ne"] + peak_alt = float(ne.alt_km[ne.argmax("alt_km")]) + assert 150 < peak_alt < 600, f"{ds_fixture} Ne peak at {peak_alt:.0f} km — outside F-region" + + +@pytest.mark.parametrize("ds_fixture", ["ds90", "ds20"]) +def test_o_plus_dominates_f_region(ds_fixture, request): + # O+ must be the dominant ion in the F-region + ds = request.getfixturevalue(ds_fixture) + f_region = ds.sel(alt_km=slice(200, 400)) + o_plus = float(f_region["O+"].mean()) + assert o_plus > float(f_region["O2+"].mean()) + assert o_plus > float(f_region["NO+"].mean()) diff --git a/tests/test_iri90.py b/tests/test_iri90.py deleted file mode 100644 index 9ac0ae9..0000000 --- a/tests/test_iri90.py +++ /dev/null @@ -1,73 +0,0 @@ -# %% -from __future__ import annotations -from datetime import datetime, UTC -import matplotlib -from glowpython2 import Iri90 -from iri20py import Iri2020 -from glowpython2.utils import alt_grid -import matplotlib.pyplot as plt -import numpy as np -# %% -usetex = False -if not usetex: - # computer modern math text - matplotlib.rcParams.update({'mathtext.fontset': 'cm'}) - -matplotlib.rc( - 'font', **{ - 'family': 'serif', - 'serif': ['Times' if usetex else 'Times New Roman'] - } -) -# for Palatino and other serif fonts use: -# rc('font',**{'family':'serif','serif':['Palatino']}) -matplotlib.rc('text', usetex=usetex) -# %% -MON = 12 -iri90 = Iri90() -iri20 = Iri2020() -_, ds90 = iri90.evaluate( - time=datetime(2015, MON, 13, 10, 0, 0, tzinfo=UTC), - lat=42.6, - lon=-71.2, - alt=alt_grid(), -) -_, ds20 = iri20.evaluate( - time=datetime(2015, MON, 13, 10, 0, 0, tzinfo=UTC), - lat=42.6, - lon=-71.2, - alt=alt_grid(), -) -# %% -ig, ax = plt.subplots(figsize=(8, 6), dpi=300) -tax = ax.twiny() -species = ['O+', 'O2+', 'N+', 'NO+'] -descs = ['O^+', 'O_2^+', 'N^+', 'NO^+'] -# species = ['N+'] -colors = ['r', 'g', 'b', 'm'] -labels = [] -lines = [] -for spec, color, desc in zip(species, colors, descs): - l21, = ds20[spec].plot(y='alt_km', ax=ax, color=color) # type: ignore - l00, = ds90[spec].plot(y='alt_km', ax=ax, linestyle='--', color=color) # type: ignore - lines.extend([l21, l00]) - labels.extend([f'IRI-20 ${desc}$', f'IRI-90 ${desc}$']) -ax.set_title('IRI-90 vs IRI-20') -l21, = ds20['Te'].plot(y='alt_km', ax=tax, color='k') # type: ignore -l00, = ds90['Te'].plot(y='alt_km', ax=tax, linestyle='--', color='k') # type: ignore -lines.extend([l21, l00]) -labels.extend(['IRI-20 $T_e$', 'IRI-90 $T_e$']) -l21, = ds20['Ti'].plot(y='alt_km', ax=tax, color='c', alpha=0.7) # type: ignore -l00, = ds90['Ti'].plot(y='alt_km', ax=tax, linestyle='--', color='c', alpha=0.7) # type: ignore -lines.extend([l21, l00]) -labels.extend(['IRI-20 $T_i$', 'IRI-90 $T_i$']) -ax.set_title('IRI-21 vs IRI-90') -ax.set_xscale('log') -ax.set_xlabel('Number Density [cm$^{-3}$]') -tax.set_xlabel('Temperature [K]') -ax.set_xlim(1e-3, None) -tax.set_xlim(100, None) -ax.legend(lines, labels, loc='upper left', fontsize='small') -plt.savefig('iri90_vs_iri20.png', dpi=300) -plt.show() -# %% diff --git a/tests/test_maxwellian.py b/tests/test_maxwellian.py new file mode 100644 index 0000000..8094af9 --- /dev/null +++ b/tests/test_maxwellian.py @@ -0,0 +1,77 @@ +import numpy as np +import pytest +import glowpython2 as glow + + +Q = 1 # erg cm-2 s-1 +Echar = 100e3 # eV + + +@pytest.fixture(scope="module") +def iono_precip(sample_time, sample_lat, sample_lon): + return glow.maxwellian(sample_time, sample_lat, sample_lon, Nbins=100, Q=Q, Echar=Echar) + + +@pytest.fixture(scope="module") +def iono_no_precip(sample_time, sample_lat, sample_lon): + return glow.no_precipitation(sample_time, sample_lat, sample_lon, Nbins=100) + + +# --- output structure --- + +def test_has_alt_coord(iono_precip): + assert "alt_km" in iono_precip.coords + + +def test_has_energy_coord(iono_precip): + assert "energy" in iono_precip.coords + + +def test_has_ver(iono_precip): + assert "ver" in iono_precip + assert "wavelength" in iono_precip["ver"].dims + + +def test_has_precip(iono_precip): + assert "precip" in iono_precip + + +# --- physical sanity --- + +def test_temperatures_positive(iono_precip): + for var in ("Te", "Ti", "Tn"): + assert np.all(iono_precip[var].values > 0) + + +def test_ver_valid_range(iono_precip): + ver = iono_precip["ver"].sel(alt_km=slice(75, None)).values + assert np.all(ver >= 0) + + +def test_ver_fill_below_boundary(iono_precip): + # Below ~72 km VER should be NaN or numerically negligible (< 1e-10 cm-3 s-1) + ver_low = iono_precip["ver"].sel(alt_km=slice(None, 72)).values + assert np.all(np.isnan(ver_low) | (np.abs(ver_low) < 1e-10)) + + +def test_precip_flux_nonzero(iono_precip): + assert np.any(iono_precip["precip"].values > 0) + + +def test_precip_increases_ver(iono_precip, iono_no_precip): + total_precip = float(iono_precip["ver"].sum()) + total_no_precip = float(iono_no_precip["ver"].sum()) + assert total_precip > total_no_precip + + +def test_auroral_4278_has_signal(iono_precip): + # N2+ first-negative (4278 A) is a direct auroral tracer — must be clearly non-zero + ver = iono_precip["ver"].sel(wavelength="4278", alt_km=slice(80, None)) + assert float(ver.max()) > 1.0 + + +def test_5577_brighter_with_precip(iono_precip, iono_no_precip): + # OI 5577 A should be at least 10x brighter with precipitation than without + v_precip = float(iono_precip["ver"].sel(wavelength="5577").max()) + v_no_precip = float(iono_no_precip["ver"].sel(wavelength="5577").max()) + assert v_precip > v_no_precip * 10 diff --git a/tests/test_msis.py b/tests/test_msis.py new file mode 100644 index 0000000..697d2b9 --- /dev/null +++ b/tests/test_msis.py @@ -0,0 +1,94 @@ +import numpy as np +import pytest +from glowpython2 import NrlMsis00 +from glowpython2.utils import alt_grid +from msis21py import NrlMsis21 + +SPECIES = ["O", "O2", "N2", "NO"] + + +@pytest.fixture(scope="module") +def ds00(sample_time, sample_lat, sample_lon): + return NrlMsis00().evaluate(time=sample_time, lat=sample_lat, lon=sample_lon, alt=alt_grid()) + + +@pytest.fixture(scope="module") +def ds21(sample_time, sample_lat, sample_lon): + return NrlMsis21().evaluate(time=sample_time, lat=sample_lat, lon=sample_lon, alt=alt_grid()) + + +# --- output structure --- + +@pytest.mark.parametrize("ds_fixture", ["ds00", "ds21"]) +def test_has_alt_coord(ds_fixture, request): + ds = request.getfixturevalue(ds_fixture) + assert "alt_km" in ds.coords + + +@pytest.mark.parametrize("ds_fixture,species", [ + (f, s) for f in ["ds00", "ds21"] for s in SPECIES +]) +def test_has_species(ds_fixture, species, request): + ds = request.getfixturevalue(ds_fixture) + assert species in ds, f"{ds_fixture} missing {species}" + + +@pytest.mark.parametrize("ds_fixture", ["ds00", "ds21"]) +def test_has_temperature(ds_fixture, request): + ds = request.getfixturevalue(ds_fixture) + assert "Tn" in ds + + +# --- physical sanity --- + +@pytest.mark.parametrize("ds_fixture", ["ds00", "ds21"]) +def test_densities_nonnegative(ds_fixture, request): + ds = request.getfixturevalue(ds_fixture) + for spec in SPECIES: + assert np.all(ds[spec].values >= 0), f"{ds_fixture} {spec} has negative values" + + +@pytest.mark.parametrize("ds_fixture", ["ds00", "ds21"]) +def test_temperature_positive(ds_fixture, request): + ds = request.getfixturevalue(ds_fixture) + assert np.all(ds["Tn"].values > 0) + + +@pytest.mark.parametrize("ds_fixture", ["ds00", "ds21"]) +def test_n2_dominates_low_alt(ds_fixture, request): + ds = request.getfixturevalue(ds_fixture) + low = ds.sel(alt_km=slice(None, 150)) + assert float(low["N2"].mean()) > float(low["O"].mean()) + + +@pytest.mark.parametrize("ds_fixture", ["ds00", "ds21"]) +def test_alt_range(ds_fixture, request): + ds = request.getfixturevalue(ds_fixture) + alt = ds.coords["alt_km"].values + assert alt.min() >= 50 + assert alt.max() <= 1100 + + +@pytest.mark.parametrize("ds_fixture", ["ds00", "ds21"]) +def test_n2_density_magnitude(ds_fixture, request): + # N2 at 100 km is well-constrained; should be ~10^12-10^14 cm^-3 + ds = request.getfixturevalue(ds_fixture) + n2 = float(ds["N2"].sel(alt_km=100, method="nearest")) + assert 1e11 < n2 < 1e14 + + +@pytest.mark.parametrize("ds_fixture", ["ds00", "ds21"]) +def test_thermosphere_warmer_than_mesosphere(ds_fixture, request): + # Thermospheric heating means Tn at 400 km > Tn at 80 km + ds = request.getfixturevalue(ds_fixture) + tn_meso = float(ds["Tn"].sel(alt_km=80, method="nearest")) + tn_thermo = float(ds["Tn"].sel(alt_km=400, method="nearest")) + assert tn_thermo > tn_meso + + +@pytest.mark.parametrize("ds_fixture", ["ds00", "ds21"]) +def test_o_dominates_above_200km(ds_fixture, request): + # Atomic O should be the dominant neutral above ~200 km + ds = request.getfixturevalue(ds_fixture) + high = ds.sel(alt_km=slice(200, None)) + assert float(high["O"].mean()) > float(high["N2"].mean()) diff --git a/tests/test_msis00.py b/tests/test_msis00.py deleted file mode 100644 index d75217c..0000000 --- a/tests/test_msis00.py +++ /dev/null @@ -1,63 +0,0 @@ -# %% -from datetime import datetime, UTC -from glowpython2 import NrlMsis00, Msis00Settings -from msis21py import NrlMsis21 -from glowpython2.utils import alt_grid -import matplotlib -import matplotlib.pyplot as plt -import numpy as np -# %% -usetex = False -if not usetex: - # computer modern math text - matplotlib.rcParams.update({'mathtext.fontset': 'cm'}) - -matplotlib.rc( - 'font', **{ - 'family': 'serif', - 'serif': ['Times' if usetex else 'Times New Roman'] - } -) -# for Palatino and other serif fonts use: -# rc('font',**{'family':'serif','serif':['Palatino']}) -matplotlib.rc('text', usetex=usetex) -# %% -MON = 1 -msis00 = NrlMsis00() -msis21 = NrlMsis21() -ds00 = msis00.evaluate( - time=datetime(2022, MON, 13, 10, 0, 0, tzinfo=UTC), - lat=0.0, - lon=0.0, - alt=alt_grid(), -) -ds21 = msis21.evaluate( - time=datetime(2022, MON, 13, 10, 0, 0, tzinfo=UTC), - lat=0.0, - lon=0.0, - alt=alt_grid(), -) -# %% -fig, ax = plt.subplots(figsize=(8, 6), dpi=300) -tax = ax.twiny() -species = ['O', 'O2', 'N2', 'NO'] -descs = ['O', 'O$_2$', 'N$_2$', 'NO'] -colors = ['r', 'g', 'b', 'm'] -labels = [] -lines = [] -for spec, color, desc in zip(species, colors, descs): - l21, = ds21[spec].plot(y='alt_km', ax=ax, color=color) - l00, = ds00[spec].plot(y='alt_km', ax=ax, linestyle='--', color=color) - lines.extend([l21, l00]) - labels.extend([f'MSIS21 {desc}', f'MSIS00 {desc}']) -ax.set_title('MSIS00 vs MSIS21') -l21, = ds21['Tn'].plot(y='alt_km', ax=tax, color='k') -l00, = ds00['Tn'].plot(y='alt_km', ax=tax, linestyle='--', color='k') -lines.extend([l21, l00]) -labels.extend(['MSIS21 $T_n$', 'MSIS00 $T_n$']) -ax.set_title('MSIS00 vs MSIS21') -ax.set_xscale('log') -ax.legend(lines, labels) -plt.savefig('msis00_vs_msis21.png', dpi=300) -plt.show() -# %% diff --git a/tests/test_no_precip.py b/tests/test_no_precip.py new file mode 100644 index 0000000..7f6f967 --- /dev/null +++ b/tests/test_no_precip.py @@ -0,0 +1,82 @@ +import numpy as np +import pytest +import glowpython2 as glow + + +@pytest.fixture(scope="module") +def iono(sample_time, sample_lat, sample_lon): + return glow.no_precipitation(sample_time, sample_lat, sample_lon, Nbins=100) + + +# --- output structure --- + +def test_has_alt_coord(iono): + assert "alt_km" in iono.coords + + +def test_alt_range(iono): + alt = iono.coords["alt_km"].values + assert alt.min() >= 50 + assert alt.max() <= 1100 + + +def test_has_neutral_species(iono): + for var in ("O", "O2", "N2", "NO"): + assert var in iono, f"Missing neutral species: {var}" + + +def test_has_ion_species(iono): + for var in ("O+", "O2+", "NO+", "NeIn"): + assert var in iono, f"Missing ion species: {var}" + + +def test_has_temperatures(iono): + for var in ("Te", "Ti", "Tn"): + assert var in iono, f"Missing temperature: {var}" + + +def test_has_ver(iono): + assert "ver" in iono + assert "wavelength" in iono["ver"].dims + + +def test_has_precip(iono): + assert "precip" in iono + + +# --- physical sanity --- + +def test_neutral_densities_positive(iono): + for var in ("O", "O2", "N2", "NO"): + assert np.all(iono[var].values >= 0), f"{var} has negative values" + + +def test_temperatures_positive(iono): + for var in ("Te", "Ti", "Tn"): + assert np.all(iono[var].values > 0), f"{var} has non-positive values" + + +def test_electron_density_positive(iono): + assert np.all(iono["NeIn"].values >= 0) + + +def test_ver_valid_range(iono): + # VER must be non-negative above 75 km + ver = iono["ver"].sel(alt_km=slice(75, None)).values + assert np.all(ver >= 0) + + +def test_ver_fill_below_boundary(iono): + # Below ~72 km GLOW outputs NaN or 0 for VER. + # Must be one of those, not arbitrary negatives + ver_low = iono["ver"].sel(alt_km=slice(None, 72)).values + assert np.all(np.isnan(ver_low) | (ver_low == 0)) + + +def test_no_precip_flux_is_zero(iono): + assert np.all(iono["precip"].values == 0) + + +def test_n2_dominates_at_low_alt(iono): + low = iono.sel(alt_km=slice(None, 150)) + assert float(low["N2"].mean()) > float(low["O"].mean())