Split __init__.py into submodules and remove the nahar handler (#72)

Splits the 1400-line `__init__.py` into focused submodules, removes the
retired nahar handler, and applies a round of cleanups. No behaviour
change: the generated ARTIS files are byte-identical to `main` for every
test set.

## Module split

`__init__.py` becomes a re-export facade (with `__all__`), so
`artisatomic.name` access from the readers, the tests, and the
`makeartisatomicfiles` entry point is unchanged. The implementation now
lives in:

| module | responsibility |
| --- | --- |
| `base.py` | element data, physical constants, shared helpers
(`log_and_print`, `xopen_check_extension`, `path_for_log`,
`parallel_map`, `ion_log_path`, …) |
| `levelnames.py` | level-name parsing: `interpret_configuration`,
`get_parity_from_config` |
| `ionhandlers.py` | which ions to process and which source reads each |
| `phixs.py` | cross-section downsampling and the hydrogenic estimate |
| `iondata.py` | `IonData` and the per-ion reading logic |
| `output.py` | the four ARTIS output-file writers |
| `cli.py` | argument parsing, `main`, the per-element loop |

## Nahar handler removed

Deletes `readnahardata.py`, the `"nahar"` branch of `read_ion_data`, the
`nahar_core_states` / `nahar_configurations` fields of `IonData`, and
the Nahar target-fraction path, plus the `atomic-data-nahar/` and unused
`recombination_rates/` folders. `makerecombratefile` keeps its optional
`.rrc` glob (now finding nothing and falling through to its other
sources).

## Cleanups

- Physical constants are defined once in `base.py`; the duplicate copies
in five readers are gone, as is a second import-time parse of
`atomic_properties.txt` in `readhillierdata`.
- `IonData` is a `dataclass(slots=True)` and the photoionisation target
fractions are filled in place, replacing `NamedTuple._replace()`.
- The Hillier-specific target fallback moved out of `write_output_files`
into `resolve_photoion_targetfractions`, so `output.py` no longer
imports any reader.
- Dead code removed: an always-true guard and a redundant `else` in the
cmfgen branch, the unused `chunks()` helper, ~30 lines of commented-out
alternatives in the downsampler.

## Verification

- 18 unit tests pass; ruff, ty, and the pre-commit hooks are clean (mypy
unchanged at its 3 pre-existing errors).
- All five `tests/*/` sets were regenerated after every commit.
`cmfgen`, `floers25`, and `jplt` match their committed checksums;
`kurucz` and `qub` produce byte-identical output to `main` when both are
run against the same local data (their committed checksums differ from
local runs both before and after this branch).
- `tests/` is unchanged from `main`, including checksums.

🤖 Generated with [Claude Code](https://claude.com/claude-code)

---------

Co-authored-by: Claude Fable 5 <noreply@anthropic.com>

Author: lukeshingles

.gitignore

@@ -14,7 +14,6 @@ output.txt
 
 *.nosync
 
-recombination_rates/phixs.pdf
 .DS_Store
 .idea
 .tags*

.sourcery.yaml

@@ -28,6 +28,8 @@ rule_settings:
         - extract-duplicate-method
         - code_clarification
         - low-code-quality
+        - no-loop-in-tests
+        - no-conditionals-in-tests
     rule_types:
         - refactoring
         - suggestion

artisatomic/__init__.py

@@ -0,0 +1,198 @@
+"""Element data, physical constants, and small utilities shared by the data-source readers."""
+
+import contextlib
+import multiprocessing as mp
+import sys
+import typing as t
+from collections.abc import Callable
+from collections.abc import Iterable
+from functools import lru_cache
+from pathlib import Path
+
+import pandas as pd
+import polars as pl
+
+PYDIR = Path(__file__).parent.resolve()
+atomicdata = pd.read_csv(PYDIR / "atomic_properties.txt", sep=r"\s+", comment="#")
+atomicdata = atomicdata.apply(lambda x: x.fillna(x.number / 0.45), axis=1)  # estimate unknown atomic mass as Z / 0.45
+elsymbols = ["n", *list(atomicdata["symbol"].values)]
+atomic_weights = ["n", *list(atomicdata["mass"].values)]
+
+roman_numerals = (
+    "",
+    "I",
+    "II",
+    "III",
+    "IV",
+    "V",
+    "VI",
+    "VII",
+    "VIII",
+    "IX",
+    "X",
+    "XI",
+    "XII",
+    "XIII",
+    "XIV",
+    "XV",
+    "XVI",
+    "XVII",
+    "XVIII",
+    "XIX",
+    "XX",
+)
+
+# the single copy of each constant for the whole package
+ryd_to_ev = 13.605693122994232
+hc_in_ev_cm = 0.0001239841984332003
+hc_in_ev_angstrom = 12398.419843320025
+h_in_ev_seconds = 4.135667696923859e-15
+
+
+def split_element_ionstage_str(ionstr: str) -> tuple[int, int]:
+    """Split a string like 'FeII' into (atomic_number, ion_stage).
+
+    Splitting on `ionstr.rstrip("IVX")` destroys the symbols of the elements whose symbols are
+    made only of those letters: V (vanadium) and I (iodine). Instead find the split point where
+    the prefix is an element symbol and the suffix is a Roman numeral. Element symbols have a
+    lowercase second letter and Roman numerals are uppercase, so the match is unambiguous.
+    """
+    for splitpos in range(1, len(ionstr)):
+        elsym, ion_stage_roman = ionstr[:splitpos], ionstr[splitpos:]
+        if elsym in elsymbols and ion_stage_roman in roman_numerals[1:]:
+            return elsymbols.index(elsym), roman_numerals.index(ion_stage_roman)
+
+    msg = f"Could not split '{ionstr}' into an element symbol and a Roman numeral ion stage"
+    raise ValueError(msg)
+
+
+# The id-keyed transition columns write_transition_data() needs, so an empty frame still carries
+# them. Name-keyed frames get lowerlevel/upperlevel from add_level_ids_forbidden() instead.
+empty_transitions_schema = pl.Schema({"lowerlevel": pl.Int64, "upperlevel": pl.Int64, "A": pl.Float64})
+
+
+def leveltuples_to_pldataframe(energy_levels) -> pl.DataFrame:
+    """Convert a list of level tuples (or a DataFrame) into a DataFrame with a zero-based levelid column.
+
+    Level ids are zero-based everywhere in memory; the 1-based numbering of the output files is
+    applied by the write_*() functions.
+    """
+    dflevels = energy_levels if isinstance(energy_levels, pl.DataFrame) else pl.DataFrame(energy_levels)
+
+    if "levelid" not in dflevels.columns:
+        dflevels = dflevels.with_row_index(name="levelid")
+
+    dflevels = dflevels.with_columns(pl.col("levelid").cast(pl.Int64))
+
+    # the frame is indexed by level id elsewhere, so a reader-supplied levelid must be contiguous
+    # and zero-based. Not an assert: input validation must survive python -O.
+    if not dflevels["levelid"].equals(pl.int_range(dflevels.height, dtype=pl.Int64, eager=True)):
+        msg = "level ids must be contiguous and start at zero"
+        raise ValueError(msg)
+
+    return dflevels
+
+
+def ion_log_path(log_folder: str | Path, atomic_number: int, ion_stage: int) -> Path:
+    """Path of the per-ion log file, written by the reading pass and appended to by the writing pass."""
+    return Path(log_folder, f"{elsymbols[atomic_number].lower()}{ion_stage:d}.txt")
+
+
+def log_and_print(flog, strout):
+    """Write a line to both stdout and this ion's log file."""
+    print(strout)
+    flog.write(strout + "\n")
+
+
+def path_for_log(filepath: str | Path) -> str:
+    """Render an input data path relative to the repository root where possible.
+
+    The log files are compared by checksum in CI, so an absolute path would make them depend on
+    where the repository happens to be checked out. Paths outside the repository (some readers
+    load data from elsewhere) are returned unchanged.
+    """
+    try:
+        return str(Path(filepath).resolve().relative_to(PYDIR.parent))
+    except ValueError:
+        return str(filepath)
+
+
+def isfloat(value: t.Any) -> bool:
+    """Whether a string parses as a float, accepting Fortran's D exponent (1.5D-3)."""
+    try:
+        float(value.replace("D", "E"))
+    except ValueError:
+        return False
+
+    return True
+
+
+compression_extensions = ("", ".zst", ".gz", ".xz")
+
+
+def find_file_check_extension(filename: str | Path) -> Path | None:
+    """Find a data file by its plain name, accepting any of the compressed variants of that name.
+
+    Returns None if neither the plain name nor any compressed form exists, so that callers which
+    treat a missing file as "no data for this ion" can say so without opening it.
+    """
+    return next((path for ext in compression_extensions if (path := Path(f"{filename}{ext}")).is_file()), None)
+
+
+def xopen_check_extension(filename: str | Path, **kwargs: t.Any) -> t.IO[t.Any]:
+    """Open a data file, trying the compressed variants of the name if it does not exist.
+
+    The data sets ship some files compressed and some not, and which ones varies between
+    downloads, so callers name the plain file and this finds whichever form is present.
+    """
+    from xopen import xopen
+
+    filepath = find_file_check_extension(filename)
+    if filepath is None:
+        filepaths = [f"{filename}{ext}" for ext in compression_extensions]
+        msg = f"Could not find any of the following files:\n  {'\n  '.join(filepaths)}."
+        raise FileNotFoundError(msg)
+
+    return xopen(filepath, **kwargs)
+
+
+@lru_cache(maxsize=1)
+def get_nist_ionization_energies_ev() -> dict[tuple[int, int], float]:
+    """Get a dictionary where dictioniz[(atomic_number, ion_sage)] = ionization_energy_ev."""
+    dfnist = pd.read_csv(
+        PYDIR / "nist_ionization.txt",
+        sep="\t",
+        usecols=["At. num", "Ion Charge", "Ionization Energy (a) (eV)"],
+    )
+
+    dictioniz = {}
+    for atomic_number, ion_charge, ioniz_ev in dfnist[
+        ["At. num", "Ion Charge", "Ionization Energy (a) (eV)"]
+    ].itertuples(index=False):
+        with contextlib.suppress(ValueError):
+            ion_stage = int(ion_charge) + 1
+            dictioniz[int(atomic_number), ion_stage] = ioniz_ev
+    return dictioniz
+
+
+def parallel_map[ResultType](
+    fn: Callable[..., ResultType],
+    *iterables: Iterable[t.Any],
+    **kwargs: t.Any,
+) -> list[ResultType]:
+    """Execute a parallel map with a progress bar using either multithreading (for free-threading python) or multiprocessing."""
+    # use a thread pool if we have no GIL (free threading)
+    use_multiprocessing = sys._is_gil_enabled()  # ruff: ignore[private-member-access]
+
+    if use_multiprocessing:
+        mp.set_start_method("spawn", force=True)
+        from tqdm.contrib.concurrent import process_map
+
+        results = process_map(fn, *iterables, **kwargs)  # type: ignore[arg-type] # zuban: ignore[no-untyped-call]
+    else:
+        from tqdm.contrib.concurrent import thread_map
+
+        results = thread_map(fn, *iterables, **kwargs)  # type: ignore[arg-type] # zuban: ignore[no-untyped-call]
+
+    assert isinstance(results, list)
+    return results

artisatomic/base.py

@@ -0,0 +1,135 @@
+#!/usr/bin/env python3
+# PYTHON_ARGCOMPLETE_OK
+"""Command-line entry point: build an ARTIS atomic database from the configured ions and handlers."""
+
+import argparse
+import glob
+import json
+import os
+import typing as t
+from collections.abc import Sequence
+from pathlib import Path
+
+import argcomplete
+
+from artisatomic import readhillierdata
+from artisatomic.iondata import read_ion_data
+from artisatomic.iondata import resolve_photoion_targetfractions
+from artisatomic.ionhandlers import get_ion_handlers
+from artisatomic.output import clear_files
+from artisatomic.output import write_compositionfile
+from artisatomic.output import write_output_files
+
+
+def main(args: argparse.Namespace | None = None, argsraw: Sequence[str] | None = None, **kwargs: t.Any) -> None:
+    """Write an ARTIS atomic database from the configured ions and handlers."""
+    if args is None:
+        parser = argparse.ArgumentParser(
+            formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+            description="Produce an ARTIS atomic database from published atomic data sets.",
+        )
+        parser.add_argument("-output_folder", action="store", default="artis_files", help="Folder for output files")
+        parser.add_argument(
+            "-output_folder_logs", action="store", default="atomic_data_logs", help="Folder for log files"
+        )
+        parser.add_argument(
+            "-nphixspoints", type=int, default=100, help="Number of cross section points to save in output"
+        )
+        parser.add_argument(
+            "-phixsnuincrement",
+            type=float,
+            default=0.03,
+            help="Fraction of nu_edge incremented for each cross section point",
+        )
+        parser.add_argument(
+            "-optimaltemperature",
+            type=int,
+            default=6000,
+            help=(
+                "(Electron and excitation) temperature at which recombination rate "
+                "should be constant when downsampling cross sections"
+            ),
+        )
+        parser.add_argument(
+            "-electrontemperature",
+            type=int,
+            default=6000,
+            help="Temperature for choosing effective collision strengths",
+        )
+        parser.add_argument(
+            "--nophixs", action="store_true", help="Don't generate cross sections and write to phixsdata_v2.txt file"
+        )
+
+        parser.add_argument(
+            "-nlevels_hydrogenic_for_unknown_phixs",
+            type=int,
+            default=100,
+            help=(
+                "Consider this many of the lowest levels of any ion whose handler supplied no"
+                " cross sections at all, and estimate a hydrogenic one for each, or 0 to disable."
+                " Negative values are rejected. Fewer tables than this can result, because a level"
+                " at or above the ionization energy is skipped but still counts towards the limit."
+                " An ion with even one cross section from its data source is left untouched, so"
+                " this never replaces or extends measured data. Excludes the top ion, which has no"
+                " upper ion to photoionise to."
+            ),
+        )
+
+        parser.set_defaults(**kwargs)
+        argcomplete.autocomplete(parser)
+        args = parser.parse_args(argsraw)
+
+    # 0 is the way to switch the estimate off, so a negative value is a typo rather than a
+    # quieter way of saying the same thing
+    if args.nlevels_hydrogenic_for_unknown_phixs < 0:
+        msg = f"-nlevels_hydrogenic_for_unknown_phixs must not be negative, got {args.nlevels_hydrogenic_for_unknown_phixs}"
+        raise ValueError(msg)
+
+    ion_handlers = get_ion_handlers()
+
+    assert len(ion_handlers) > 0
+    readhillierdata.read_hyd_phixsdata()
+
+    os.makedirs(args.output_folder, exist_ok=True)
+
+    log_folder = Path(args.output_folder) / args.output_folder_logs
+    if log_folder.exists():
+        # delete any existing log files
+        logfiles = glob.glob(os.path.join(log_folder, "*.txt"))
+        for logfile in logfiles:
+            Path(logfile).unlink(missing_ok=True)
+            print("deleting", logfile)
+    else:
+        os.makedirs(log_folder, exist_ok=True)
+
+    with Path(log_folder, "artisatomicionhandlers.json").open("w", encoding="utf-8") as f:
+        json.dump(obj=ion_handlers, fp=f)
+    write_compositionfile(ion_handlers, args)
+    clear_files(args)
+    process_files(ion_handlers, args)
+
+
+def process_files(ion_handlers: list[tuple[int, list[tuple[int, str]]]], args: argparse.Namespace) -> None:
+    """Read every configured ion and append it to the output files, one element at a time.
+
+    Ion stages are processed in ascending order so that each ion's photoionisation targets, which
+    are levels of the next ion up, are already known, and so the top ion can be identified and
+    given no cross sections.
+    """
+    for atomic_number, listions in ion_handlers:
+        if not listions:
+            continue
+
+        iondatalist = [
+            read_ion_data(atomic_number, ion_stage_entry, is_top_ion=(i == len(listions) - 1), args=args)
+            for i, ion_stage_entry in enumerate(listions)
+        ]
+
+        if not args.nophixs:
+            resolve_photoion_targetfractions(iondatalist)
+
+        write_output_files(atomic_number, iondatalist, args)
+
+
+if __name__ == "__main__":
+    main()

artisatomic/cli.py

@@ -9,8 +9,6 @@
 
 import artisatomic
 
-hc_in_ev_cm = 0.0001239841984332003
-
 
 class EnergyLevel(t.NamedTuple):
     """A ground state read from the NIST table."""