diff --git a/CHANGELOG.md b/CHANGELOG.md
index ffd25625..d1fdc054 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -223,6 +223,51 @@ returned.
   disable/re-enable cycle and are re-attached when magnetism is enabled
   again. `update_layer` also accepts the magnetism keys one at a time.
 
+## ORSO file handling
+
+- Binary ORSO (`.orb`, NeXus/HDF5) files are read and written alongside
+  `.ort` text files. The format is detected from the ORSO banner line or
+  the HDF5 magic bytes, not the file extension, so a file that carries
+  the banner but fails to parse now raises instead of being silently
+  re-read as plain text (which dropped the whole header, polarization
+  included). `.orb` support needs `h5py`, available as the new `orb`
+  extra; `orsopy` is pinned to `>=1.2`.
+- New ORSO export. `Project.save_experiment_as_orso(path, index=None)`
+  writes an experiment (`.ort` or `.orb`), with the model, when set,
+  serialized as `data_source.sample.model`. Backed by the new
+  `save_orso_experiment`, `orso_datasets_from_experiment` and
+  `sample_to_orso_model` in `easyreflectometry.orso_utils`. A polarized
+  experiment becomes one file with one `data_set:` block per spin
+  channel. `Model.as_orso` now returns the ORSO model-language
+  dictionary (slab representation) rather than the internal `as_dict`.
+- Repeating multilayers survive a round trip. Loading resolves the ORSO
+  stack with `resolve_stack()` instead of flattening it, so a sub-stack
+  keeps its repetition count and comes back as a `RepeatingMultilayer`;
+  export writes it with the inline `N ( ... )` stack syntax.
+- Units declared in the file are honoured: `Qz` in `1/nm`, lengths in
+  `nm` (the ORSO default) and SLDs in `1/nm^2` are converted on load,
+  instead of being read as angstrom-based numbers.
+- Resolution and error columns are read more carefully. A column
+  declared `value_is: FWHM` is converted to sigma on load, `nan` entries
+  in `sQz` are filled by interpolating over the valid points, and
+  partially missing error columns warn rather than propagating `nan`
+  into a fit. Stored `Pointwise` resolutions remain variances, so saved
+  projects round-trip without migration.
+- `Project.load_polarized_experiment_from_file(path)` loads a polarized
+  experiment from a single multi-dataset ORSO file, classifying each
+  `data_set:` block by its own `instrument_settings.polarization`
+  header. Only `pp/pm/mp/mm` are mapped; a file with an unmappable or
+  duplicated channel raises rather than guessing. Supported by the new
+  `channel_from_orso_polarization` and
+  `detect_polarization_channels_per_dataset` in
+  `easyreflectometry.data`.
+- New `easyreflectometry.data.dataset_from_datagroup` builds a
+  `DataSet1D` from one dataset of an already-loaded `DataGroup`, and
+  keeps the parsed ORSO header on the dataset as `orso_header` so
+  exporters can reuse the original provenance. `load_as_dataset` and the
+  project loaders accept a pre-loaded `DataGroup`, so importing a file
+  no longer parses it three or four times.
+
 ## Documentation
 
 - The documentation is now MkDocs (Material) only. The legacy Sphinx
diff --git a/pixi.lock b/pixi.lock
index e457f388..b2006eae 100644
--- a/pixi.lock
+++ b/pixi.lock
@@ -8316,7 +8316,7 @@ packages:
   - asteval
   - bumps
   - easyscience @ git+https://github.com/easyscience/easyscience.git@sampler-engine-structure-280
-  - orsopy
+  - orsopy>=1.2
   - plotly
   - pooch
   - refl1d>=1.0.0
@@ -8358,6 +8358,7 @@ packages:
   - spdx-headers ; extra == 'dev'
   - validate-pyproject[all] ; extra == 'dev'
   - versioningit ; extra == 'dev'
+  - h5py ; extra == 'orb'
   requires_python: '>=3.11'
 - pypi: git+https://github.com/easyscience/easyscience.git?rev=sampler-engine-structure-280#3e598107beb2de6820e0e4c0a68da8e3c9af4ae4
   name: easyscience
diff --git a/pyproject.toml b/pyproject.toml
index f646fd75..eaf383ea 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -27,7 +27,8 @@ dependencies = [
   'scipp',
   'refnx',
   'refl1d>=1.0.0',
-  'orsopy',
+  # >=1.2: model language resolve_to_blocks and nexus (.orb) I/O verified in 1.2.2
+  'orsopy>=1.2',
   'svglib<1.6 ; platform_system=="Linux" or sys_platform == "darwin"',
   'xhtml2pdf',
   'bumps',
@@ -37,6 +38,8 @@ dependencies = [
 ]
 
 [project.optional-dependencies]
+# Binary ORSO (.orb / NeXus) read and write support (orsopy uses h5py for it)
+orb = ['h5py']
 dev = [
   'GitPython',                       # Interact with Git repositories
   'build',                           # Building the package
diff --git a/src/easyreflectometry/data/__init__.py b/src/easyreflectometry/data/__init__.py
index e63470eb..4650b63d 100644
--- a/src/easyreflectometry/data/__init__.py
+++ b/src/easyreflectometry/data/__init__.py
@@ -3,18 +3,24 @@
 
 from .data_store import DataSet1D
 from .data_store import ProjectData
+from .measurement import dataset_from_datagroup
 from .measurement import load
 from .measurement import load_as_dataset
 from .measurement import merge_datagroups
 from .polarized import PolarizedDataSet
+from .polarized import channel_from_orso_polarization
 from .polarized import detect_polarization_channel
+from .polarized import detect_polarization_channels_per_dataset
 
 __all__ = [
     'load',
     'load_as_dataset',
+    'dataset_from_datagroup',
     'merge_datagroups',
     'ProjectData',
     'DataSet1D',
     'PolarizedDataSet',
+    'channel_from_orso_polarization',
     'detect_polarization_channel',
+    'detect_polarization_channels_per_dataset',
 ]
diff --git a/src/easyreflectometry/data/measurement.py b/src/easyreflectometry/data/measurement.py
index 12117162..bd3c7432 100644
--- a/src/easyreflectometry/data/measurement.py
+++ b/src/easyreflectometry/data/measurement.py
@@ -3,6 +3,7 @@
 
 
 import os
+from typing import Optional
 from typing import TextIO
 from typing import Union
 
@@ -10,40 +11,97 @@
 import scipp as sc
 
 from easyreflectometry.data import DataSet1D
+from easyreflectometry.orso_utils import is_orso_file
 from easyreflectometry.orso_utils import load_data_from_orso_file
 
 
 def load(fname: Union[TextIO, str]) -> sc.DataGroup:
-    """Load data from an ORSO .ort file.
+    """Load data from an ORSO file (.ort/.orb) or a plain text file.
+
+    The discriminator is the ORSO banner line (or the HDF5 magic for binary
+    files), **not** the file extension: a file carrying the banner that fails
+    to parse raises instead of being silently re-read as plain text (which
+    would drop the entire header, including polarization).
 
     Parameters
     ----------
     fname : Union[TextIO, str]
         The file to be read.
+
+    Returns
+    -------
+    sc.DataGroup
+        The loaded data.
     """
-    try:
+    if is_orso_file(str(fname)):
         return load_data_from_orso_file(fname)
-    except (IndexError, ValueError):
-        return _load_txt(fname)
+    return _load_txt(fname)
 
 
-def load_as_dataset(fname: Union[TextIO, str]) -> DataSet1D:
-    """Load data from an ORSO .ort file as a DataSet1D."""
-    data_group = load(fname)
-    basename = os.path.splitext(os.path.basename(fname))[0]
-    data_name = 'R_' + basename
-    coords_name = 'Qz_' + basename
-    coords_name = list(data_group['coords'].keys())[0] if coords_name not in data_group['coords'] else coords_name
-    data_name = list(data_group['data'].keys())[0] if data_name not in data_group['data'] else data_name
+def dataset_from_datagroup(data_group: sc.DataGroup, data_key: Optional[str] = None) -> DataSet1D:
+    """Build a DataSet1D from one dataset of a loaded DataGroup.
+
+    The ORSO header (when present) is attached to the returned dataset as the
+    ``orso_header`` attribute (a plain dict), so exporters can reuse the
+    original ``data_source``/``reduction`` provenance.
+
+    Parameters
+    ----------
+    data_group : sc.DataGroup
+        A DataGroup as returned by :func:`load`.
+    data_key : Optional[str], optional
+        The data entry to use (e.g. ``'R_0'``). By default, the first entry.
+
+    Returns
+    -------
+    DataSet1D
+        The dataset.
+    """
+    if data_key is None:
+        data_key = list(data_group['data'].keys())[0]
+    coords_key = 'Qz_' + data_key[len('R_') :]
+    if coords_key not in data_group['coords']:
+        coords_key = list(data_group['coords'].keys())[0]
     dataset = DataSet1D(
-        x=data_group['coords'][coords_name].values,
-        y=data_group['data'][data_name].values,
-        ye=data_group['data'][data_name].variances,
-        xe=data_group['coords'][coords_name].variances,
+        x=data_group['coords'][coords_key].values,
+        y=data_group['data'][data_key].values,
+        ye=data_group['data'][data_key].variances,
+        xe=data_group['coords'][coords_key].variances,
     )
+    header = None
+    if 'attrs' in data_group and data_key in data_group['attrs']:
+        try:
+            header = data_group['attrs'][data_key]['orso_header'].values
+        except (KeyError, AttributeError):
+            header = None
+    dataset.orso_header = header
     return dataset
 
 
+def load_as_dataset(fname: Union[TextIO, str], data_group: Optional[sc.DataGroup] = None) -> DataSet1D:
+    """Load data from an ORSO .ort file as a DataSet1D.
+
+    Parameters
+    ----------
+    fname : Union[TextIO, str]
+        The file to be read.
+    data_group : Optional[sc.DataGroup], optional
+        Pre-loaded DataGroup for *fname* (avoids re-parsing the file).
+        By default, None.
+
+    Returns
+    -------
+    DataSet1D
+        The (first) dataset in the file.
+    """
+    if data_group is None:
+        data_group = load(fname)
+    basename = os.path.splitext(os.path.basename(fname))[0]
+    data_name = 'R_' + basename
+    data_name = list(data_group['data'].keys())[0] if data_name not in data_group['data'] else data_name
+    return dataset_from_datagroup(data_group, data_key=data_name)
+
+
 def extract_orso_title(data_group: sc.DataGroup, data_name: str) -> str | None:
     """Extract orso title."""
     try:
diff --git a/src/easyreflectometry/data/polarized.py b/src/easyreflectometry/data/polarized.py
index ae8fcac8..75eaa39c 100644
--- a/src/easyreflectometry/data/polarized.py
+++ b/src/easyreflectometry/data/polarized.py
@@ -244,6 +244,60 @@ def detect_polarization_channel(path: str) -> Optional[PolarizationChannel]:
     return _channel_from_filename(path)
 
 
+def channel_from_orso_polarization(polarization) -> Optional[PolarizationChannel]:
+    """Map an ORSO ``instrument_settings.polarization`` value to a spin channel.
+
+    Parameters
+    ----------
+    polarization :
+        The header value (orsopy ``Polarization`` enum, string, or None).
+
+    Returns
+    -------
+    Optional[PolarizationChannel]
+        The mapped channel, or None for absent/unmapped values (``po``, ``mo``,
+        ``op``, ``om``, ``unpolarized``, ``vector`` are deliberately unmapped).
+    """
+    if polarization is None:
+        return None
+    value = getattr(polarization, 'value', polarization)
+    return _ORSO_POLARIZATION_TO_CHANNEL.get(str(value).lower())
+
+
+def _dataset_polarization(orso_dataset):
+    """The declared polarization of one parsed ORSO dataset, or None."""
+    try:
+        return orso_dataset.info.data_source.measurement.instrument_settings.polarization
+    except AttributeError:
+        return None
+
+
+def detect_polarization_channels_per_dataset(
+    orso_data,
+) -> list[tuple[bool, Optional[PolarizationChannel]]]:
+    """Classify every dataset of a parsed ORSO file by its own header.
+
+    Unlike :func:`detect_polarization_channel`, which reads only the first
+    dataset, this honours per-dataset ``polarization:`` overrides in
+    multi-dataset files.
+
+    Parameters
+    ----------
+    orso_data : list
+        Parsed ORSO dataset list (as returned by ``orso.load_orso``).
+
+    Returns
+    -------
+    list[tuple[bool, Optional[PolarizationChannel]]]
+        Per dataset: (header declares a polarization, mapped channel or None).
+    """
+    result = []
+    for orso_dataset in orso_data:
+        polarization = _dataset_polarization(orso_dataset)
+        result.append((polarization is not None, channel_from_orso_polarization(polarization)))
+    return result
+
+
 def _channel_from_orso_header(path: str) -> tuple[bool, Optional[PolarizationChannel]]:
     """Read the polarization of the first dataset in an ORSO file.
 
@@ -254,16 +308,15 @@ def _channel_from_orso_header(path: str) -> tuple[bool, Optional[PolarizationCha
         False when the file is unreadable or carries no polarization field.
     """
     try:
-        from orsopy.fileio import orso
+        from easyreflectometry.orso_utils import _load_orso_any
 
-        orso_data = orso.load_orso(str(path))
+        orso_data = _load_orso_any(str(path))
         polarization = orso_data[0].info.data_source.measurement.instrument_settings.polarization
     except Exception:
         return False, None
     if polarization is None:
         return False, None
-    value = getattr(polarization, 'value', polarization)
-    return True, _ORSO_POLARIZATION_TO_CHANNEL.get(str(value).lower())
+    return True, channel_from_orso_polarization(polarization)
 
 
 def _channel_from_filename(path: str) -> Optional[PolarizationChannel]:
diff --git a/src/easyreflectometry/model/model.py b/src/easyreflectometry/model/model.py
index 4cb47142..22b35f3f 100644
--- a/src/easyreflectometry/model/model.py
+++ b/src/easyreflectometry/model/model.py
@@ -344,8 +344,23 @@ def as_dict(self, skip: Optional[list[str]] = None) -> dict:
         return self.to_dict(skip=skip)
 
     def as_orso(self) -> dict:
-        """Convert the model to a dictionary suitable for ORSO."""
-        return self.as_dict()
+        """The sample as an ORSO simple-model (``sample.model``) dictionary.
+
+        Slab representation: lengths in angstrom, SLDs in 1/angstrom^2,
+        repeating multilayers via the inline ``N ( ... )`` stack syntax.
+
+        Returns
+        -------
+        dict
+            The ORSO model-language dictionary (the content of an .ort file's
+            ``data_source.sample.model`` section).
+        """
+        # Circular import if hoisted to module-top.
+        from orsopy.fileio import Header
+
+        from easyreflectometry.orso_utils import sample_to_orso_model
+
+        return Header.asdict(sample_to_orso_model(self.sample))
 
     @classmethod
     def from_dict(cls, passed_dict: dict) -> Model:
diff --git a/src/easyreflectometry/model/resolution_functions.py b/src/easyreflectometry/model/resolution_functions.py
index 9579ad66..646556f6 100644
--- a/src/easyreflectometry/model/resolution_functions.py
+++ b/src/easyreflectometry/model/resolution_functions.py
@@ -119,6 +119,12 @@ class Pointwise(ResolutionFunction):
     widths from the ``[Qz, R, sQz]`` triple loaded from a data file; the returned
     widths are consumed by the calculators (refnx ``x_err`` / refl1d ``dq``),
     which perform the actual convolution against the model.
+
+    Serialization contract: ``as_dict``/``from_dict`` store ``sQz_data_points``
+    as **variances** (sigma squared). This is deliberately unchanged by the
+    ORSO ``value_is: FWHM`` support — FWHM columns are converted to sigma at
+    load time, so stored values are always sigma squared and saved projects
+    round-trip without migration.
     """
 
     def __init__(self, q_data_points: List[np.ndarray]):
diff --git a/src/easyreflectometry/orso_utils.py b/src/easyreflectometry/orso_utils.py
index 15cb717e..75984abc 100644
--- a/src/easyreflectometry/orso_utils.py
+++ b/src/easyreflectometry/orso_utils.py
@@ -1,19 +1,45 @@
 # SPDX-FileCopyrightText: 2025 EasyScience contributors <https://github.com/easyscience>
 # SPDX-License-Identifier: BSD-3-Clause
 
+"""ORSO file support: reading and writing `.ort` (text) and `.orb` (binary) files.
+
+Reading is built on ``orsopy.fileio``.  Data columns follow the ORSO
+specification: the first four columns are Qz, R, sR, sQz (fixed order); sR and
+sQz default to sigma but a ``value_is: FWHM`` declaration is honoured and
+converted to sigma on load.  Q values declared in ``1/nm`` are converted to the
+internal ``1/angstrom`` unit.  Errors are stored as **variances** on the scipp
+arrays (sigma squared), which is also the convention used by
+:class:`~easyreflectometry.model.resolution_functions.Pointwise` and by project
+serialization (``sQz_data_points``); any exporter must convert back to sigma.
+
+The ORSO "simple model" language is resolved with units honoured (the language
+default length unit is **nm**; internal thicknesses/roughnesses are angstrom)
+and with sub-stack repetitions mapped to
+:class:`~easyreflectometry.sample.assemblies.repeating_multilayer.RepeatingMultilayer`.
+"""
+
 import logging
 import warnings
+from typing import List
+from typing import Optional
+from typing import Union
 
 import numpy as np
 import scipp as sc
 from orsopy.fileio import Header
 from orsopy.fileio import model_language
 from orsopy.fileio import orso
+from orsopy.fileio.base import Column
 from orsopy.fileio.base import ComplexValue
+from orsopy.fileio.base import ErrorColumn
+from orsopy.fileio.base import Value
+from orsopy.fileio.orso import Orso
+from orsopy.fileio.orso import OrsoDataset
 
 from easyreflectometry.data import DataSet1D
 
 from .sample.assemblies.multilayer import Multilayer
+from .sample.assemblies.repeating_multilayer import RepeatingMultilayer
 from .sample.collections.sample import Sample
 from .sample.elements.layers.layer import Layer
 from .sample.elements.materials.material import Material
@@ -22,9 +48,115 @@
 # Set up logging
 logger = logging.getLogger(__name__)
 
+# Mirrors resolution_functions.SIGMA_TO_FWHM; kept local to avoid importing the
+# model package from this low-level module (data <-> model import cycle).
+SIGMA_TO_FWHM = 2 * np.sqrt(2 * np.log(2))
+
+# The mandatory first line of a text ORSO file; the discriminator for the
+# ORSO-vs-plain-text decision (the file extension is not reliable).
+ORSO_BANNER = 'ORSO reflectivity data file'
+
+# Magic bytes of an HDF5 container -- how a binary ORSO (.orb, NeXus) file starts.
+_HDF5_MAGIC = b'\x89HDF\r\n\x1a\n'
+
+# Length unit -> angstrom conversion factors for model-language values.
+_LENGTH_UNIT_TO_ANGSTROM = {
+    None: 1.0,
+    'angstrom': 1.0,
+    'A': 1.0,
+    'nm': 10.0,
+    'um': 1.0e4,
+    'mm': 1.0e7,
+}
+
+# SLD unit -> 1/angstrom^2 conversion factors.
+_SLD_UNIT_TO_INV_ANGSTROM_SQUARED = {
+    None: 1.0,
+    '1/angstrom^2': 1.0,
+    '1/nm^2': 1.0e-2,
+}
+
+# Qz column unit -> 1/angstrom conversion factors.
+_Q_UNIT_TO_INV_ANGSTROM = {
+    None: 1.0,
+    '1/angstrom': 1.0,
+    '1/nm': 0.1,
+}
+
+# Mass density unit -> g/cm^3 conversion factors.
+_MASS_DENSITY_UNIT_TO_G_CM3 = {
+    None: 1.0,
+    'g/cm^3': 1.0,
+    'g/ml': 1.0,
+    'kg/m^3': 1.0e-3,
+}
+
+
+def is_orso_file(fname: str) -> bool:
+    """Whether *fname* is an ORSO file: text with the ORSO banner line, or HDF5 (.orb).
+
+    Parameters
+    ----------
+    fname : str
+        Path to the file.
+
+    Returns
+    -------
+    bool
+        True when the file starts with the ORSO banner or the HDF5 magic bytes.
+    """
+    try:
+        with open(fname, 'rb') as f:
+            head = f.read(128)
+    except OSError:
+        return False
+    if head.startswith(_HDF5_MAGIC):
+        return True
+    try:
+        first_line = head.decode('utf-8', errors='replace').splitlines()[0]
+    except IndexError:
+        return False
+    return first_line.lstrip().startswith('#') and ORSO_BANNER in first_line
+
+
+def _is_binary_orso(fname: str) -> bool:
+    """Whether *fname* is a binary (HDF5 / .orb) ORSO file."""
+    try:
+        with open(fname, 'rb') as f:
+            return f.read(8).startswith(_HDF5_MAGIC)
+    except OSError:
+        return False
+
+
+def _load_orso_any(fname: str) -> List[OrsoDataset]:
+    """Parse an ORSO file, text (`.ort`) or binary (`.orb`), into OrsoDataset objects.
+
+    Parameters
+    ----------
+    fname : str
+        Path to the file.
+
+    Returns
+    -------
+    List[OrsoDataset]
+        The parsed datasets.
+
+    Raises
+    ------
+    ValueError :
+        If the file cannot be parsed as ORSO (the original error is chained).
+    """
+    fname = str(fname)
+    try:
+        if _is_binary_orso(fname):
+            return orso.load_nexus(fname)
+        return orso.load_orso(fname)
+    except Exception as e:
+        raise ValueError(f'Error loading ORSO file {fname!r}: {e}') from e
+
 
 def LoadOrso(orso_data):
-    """Load a model from an ORSO file."""
+    """Load a model and data from an ORSO file (path or pre-parsed datasets)."""
 
     orso_obj = _coerce_orso_object(orso_data)
     sample = load_orso_model(orso_obj)
@@ -39,28 +171,252 @@ def _coerce_orso_object(orso_input):
             return orso_input
     except (TypeError, IndexError):
         pass
-    return orso.load_orso(orso_input)
+    return _load_orso_any(orso_input)
 
 
 def load_data_from_orso_file(fname: str) -> sc.DataGroup:
-    """Load data from an ORSO file."""
-    try:
-        orso_data = orso.load_orso(fname)
-    except Exception as e:
-        raise ValueError(f'Error loading ORSO file: {e}')
+    """Load data from an ORSO file (`.ort` text or `.orb` binary).
+
+    Parameters
+    ----------
+    fname : str
+        Path to the file.
+
+    Returns
+    -------
+    sc.DataGroup
+        A scipp DataGroup with data, coords, and attrs.
+
+    Raises
+    ------
+    ValueError :
+        If the file cannot be parsed as ORSO. Parse failures are **not**
+        swallowed here; falling back to plain-text loading is the caller's
+        decision and only valid when the file carries no ORSO banner.
+    """
+    orso_data = _load_orso_any(fname)
     return load_orso_data(orso_data)
 
 
+def _orso_dataset_key(o, index: int) -> Union[str, int]:
+    """The name a dataset is stored under: its ``data_set`` label or its index."""
+    if o.info.data_set is not None:
+        return o.info.data_set
+    return index
+
+
+def _validate_columns(columns, dataset_label) -> None:
+    """Warn when the leading four columns do not follow the ORSO layout.
+
+    The spec fixes the order Qz, R, sR, sQz; data is read by position, so a
+    file with different columns is very likely misread. The names in
+    ``info.columns`` are used as validation only.
+    """
+    expected = ('Qz', 'R', 'R', 'Qz')
+    for position, expected_name in enumerate(expected):
+        if position >= len(columns):
+            warnings.warn(
+                f'ORSO dataset {dataset_label!r} declares only {len(columns)} columns; the specification '
+                f'requires Qz, R, sR, sQz (nan-filled when unknown). Missing error columns are treated as absent.',
+                UserWarning,
+                stacklevel=3,
+            )
+            return
+        column = columns[position]
+        name = getattr(column, 'error_of', None) or getattr(column, 'name', None)
+        # Legacy (0.1 standard) files use plain columns named 'sR'/'sQz'.
+        if position >= 2 and isinstance(name, str) and name.startswith('s'):
+            name = name[1:]
+        if name != expected_name:
+            warnings.warn(
+                f'ORSO dataset {dataset_label!r} columns do not follow the specified order '
+                f'(Qz, R, sR, sQz): column {position} is {name!r}. Data is read by position '
+                f'and may be misinterpreted.',
+                UserWarning,
+                stacklevel=3,
+            )
+            return
+
+
+def _q_unit_scale(column, dataset_label) -> float:
+    """Conversion factor from the Qz column unit to the internal 1/angstrom."""
+    unit = getattr(column, 'unit', None)
+    try:
+        return _Q_UNIT_TO_INV_ANGSTROM[unit]
+    except KeyError:
+        warnings.warn(
+            f'ORSO dataset {dataset_label!r} declares Qz unit {unit!r}; expected 1/angstrom or 1/nm. '
+            f'Values are used as-is (assumed 1/angstrom).',
+            UserWarning,
+            stacklevel=3,
+        )
+        return 1.0
+
+
+def _error_column_sigma(o, position: int) -> Optional[np.ndarray]:
+    """Extract an error column as sigma values, or None when the column is absent.
+
+    Honours ``value_is: FWHM`` by converting to sigma (dividing by 2.3548...).
+    """
+    if o.data.ndim != 2 or o.data.shape[1] <= position:
+        return None
+    values = np.asarray(o.data[:, position], dtype=float)
+    if position < len(o.info.columns):
+        value_is = getattr(o.info.columns[position], 'value_is', None)
+        if value_is == 'FWHM':
+            values = values / SIGMA_TO_FWHM
+    return values
+
+
+def _clean_sqz(sqz: Optional[np.ndarray], dataset_label) -> Optional[np.ndarray]:
+    """Apply the nan policy to the sQz column.
+
+    All-nan (spec-valid "resolution unknown") returns None so that no q-variance
+    is stored and the caller falls back to the default percentage smearing.
+    Partial nan (common at the low/high-Q extremes of real reductions) is filled
+    by interpolating sigma over the valid points -- otherwise a single nan would
+    propagate through ``np.interp`` in ``Pointwise`` across the whole range.
+    """
+    if sqz is None:
+        return None
+    nan_mask = np.isnan(sqz)
+    if not nan_mask.any():
+        return sqz
+    if nan_mask.all():
+        return None
+    warnings.warn(
+        f'ORSO dataset {dataset_label!r}: {nan_mask.sum()} of {sqz.size} sQz values are nan; '
+        f'they are filled by interpolating the resolution over the valid points.',
+        UserWarning,
+        stacklevel=3,
+    )
+    valid = ~nan_mask
+    indices = np.arange(sqz.size)
+    filled = sqz.copy()
+    filled[nan_mask] = np.interp(indices[nan_mask], indices[valid], sqz[valid])
+    return filled
+
+
+def _clean_sr(sr: Optional[np.ndarray], dataset_label) -> Optional[np.ndarray]:
+    """Apply the nan policy to the sR column.
+
+    All-nan returns None (uncertainty unknown). Partial nan is kept as-is --
+    interpolating measurement uncertainties would fabricate fit weights -- but
+    is warned about, since nan weights degrade fitting.
+    """
+    if sr is None:
+        return None
+    nan_mask = np.isnan(sr)
+    if nan_mask.all() and sr.size:
+        return None
+    if nan_mask.any():
+        warnings.warn(
+            f'ORSO dataset {dataset_label!r}: {nan_mask.sum()} of {sr.size} sR values are nan; '
+            f'these points carry no uncertainty and will degrade fit weighting.',
+            UserWarning,
+            stacklevel=3,
+        )
+    return sr
+
+
+def load_orso_data(orso_data) -> sc.DataGroup:
+    """Convert parsed ORSO dataset objects into a scipp DataGroup.
+
+    Q values are converted to 1/angstrom; sR/sQz are converted to sigma when
+    declared as FWHM and stored as variances; nan-filled error columns follow
+    the policy documented on the cleaning helpers.
+
+    Parameters
+    ----------
+    orso_data : list
+        Parsed ORSO dataset list (as returned by ``orso.load_orso``).
+
+    Returns
+    -------
+    sc.DataGroup
+        A scipp DataGroup with data, coords, and attrs.
+    """
+    data = {}
+    coords = {}
+    attrs = {}
+    for i, o in enumerate(orso_data):
+        name = _orso_dataset_key(o, i)
+        _validate_columns(o.info.columns, name)
+
+        q_scale = _q_unit_scale(o.info.columns[0], name)
+        qz = np.asarray(o.data[:, 0], dtype=float) * q_scale
+        reflectivity = np.asarray(o.data[:, 1], dtype=float)
+
+        sr = _clean_sr(_error_column_sigma(o, 2), name)
+        sqz = _clean_sqz(_error_column_sigma(o, 3), name)
+        if sqz is not None:
+            sqz = sqz * q_scale
+
+        dims = [f'{o.info.columns[0].name}_{name}']
+        coords[f'Qz_{name}'] = sc.array(
+            dims=dims,
+            values=qz,
+            variances=np.square(sqz) if sqz is not None else None,
+            unit=sc.Unit('1/angstrom'),
+        )
+        r_unit = getattr(o.info.columns[1], 'unit', None)
+        try:
+            data[f'R_{name}'] = sc.array(
+                dims=dims,
+                values=reflectivity,
+                variances=np.square(sr) if sr is not None else None,
+                unit=sc.Unit(r_unit) if r_unit is not None else None,
+            )
+        except TypeError:
+            data[f'R_{name}'] = sc.array(
+                dims=dims,
+                values=reflectivity,
+                variances=np.square(sr) if sr is not None else None,
+            )
+        attrs[f'R_{name}'] = {'orso_header': sc.scalar(Header.asdict(o.info))}
+    data_group = sc.DataGroup(data=data, coords=coords, attrs=attrs)
+    return data_group
+
+
+# ---------------------------------------------------------------------------
+# Model language -> Sample
+# ---------------------------------------------------------------------------
+
+
+def _length_to_angstrom(value) -> float:
+    """Convert a model-language length (Value with unit, or bare number) to angstrom."""
+    if value is None:
+        return 0.0
+    magnitude = getattr(value, 'magnitude', value)
+    if magnitude is None:
+        return 0.0
+    unit = getattr(value, 'unit', None)
+    try:
+        factor = _LENGTH_UNIT_TO_ANGSTROM[unit]
+    except KeyError:
+        warnings.warn(
+            f'Unknown ORSO length unit {unit!r}; value used as angstrom.',
+            UserWarning,
+            stacklevel=4,
+        )
+        factor = 1.0
+    return float(magnitude) * factor
+
+
 def load_orso_model(orso_data) -> Sample:
     """Load a model from an ORSO file and return a Sample object.
 
-    The ORSO file .ort contains information about the sample, saved
-    as a simple "stack" string, e.g. 'air | m1 | SiO2 | Si'.
-    This gets parsed by the ORSO library and converted into an ORSO Dataset object.
+    The **original** ``sample.model`` object is resolved (keeping ``globals``,
+    ``materials``, ``sub_stacks`` and ``composits``), so declared units --
+    including the model-language default length unit of **nm** -- are honoured
+    and converted to the internal angstrom.  Sub-stacks with repetitions map to
+    :class:`RepeatingMultilayer`; density-defined materials stay density-defined
+    (:class:`MaterialDensity`) instead of being flattened to a numeric SLD.
 
     The stack is converted to a proper Sample structure:
     - First layer -> Superphase assembly (thickness=0, roughness=0, both fixed)
-    - Middle layers -> 'Loaded layer' Multilayer assembly (parameters enabled)
+    - Middle layers -> 'Loaded layer' Multilayer assembly (parameters enabled),
+      with repeated sub-stacks as RepeatingMultilayer assemblies
     - Last layer -> Subphase assembly (thickness=0 fixed, roughness enabled)
 
     Parameters
@@ -78,7 +434,6 @@ def load_orso_model(orso_data) -> Sample:
     Sample
         An EasyReflectometry Sample object.
     """
-    # Extract stack string and layer definitions from ORSO sample model
     sample_model = orso_data[0].info.data_source.sample.model
     if sample_model is None:
         warnings.warn(
@@ -105,29 +460,39 @@ def load_orso_model(orso_data) -> Sample:
             globals=getattr(sample_model, 'globals', None),
         )
 
-    # Try to resolve layers using different methods
-    try:
-        orso_layers = orso_sample.resolve_to_layers()
-    except ValueError:
-        orso_layers = orso_sample.resolve_stack()
+    # Resolve the original model (globals/materials/sub_stacks intact) at the
+    # stack level: resolve_stack() keeps SubStack objects (and with them the
+    # repetition counts) that resolve_to_layers()/resolve_to_blocks() flatten.
+    orso_blocks = orso_sample.resolve_stack()
 
     # Handle case where layers are not resolved correctly
-    if not orso_layers:
+    if not orso_blocks:
         raise ValueError('Could not resolve ORSO layers.')
 
-    if len(orso_layers) < 2:
+    # Plain layers still need their material resolved (what
+    # resolve_to_layers() would have done); a failure is warned about instead
+    # of silently re-resolving the whole stack a different way.
+    for block in orso_blocks:
+        if isinstance(block, model_language.Layer):
+            _generate_layer_material(block)
+
+    # The ambient (first) and substrate (last) entries must be plain layers;
+    # flatten pathological edge sub-stacks.
+    if not isinstance(orso_blocks[0], model_language.Layer):
+        warnings.warn('First ORSO stack item is a sub-stack; its layers are used directly.', UserWarning, stacklevel=2)
+        orso_blocks = list(orso_blocks[0].resolve_to_layers()) + orso_blocks[1:]
+    if not isinstance(orso_blocks[-1], model_language.Layer):
+        warnings.warn('Last ORSO stack item is a sub-stack; its layers are used directly.', UserWarning, stacklevel=2)
+        orso_blocks = orso_blocks[:-1] + list(orso_blocks[-1].resolve_to_layers())
+
+    total_layers = sum(1 if isinstance(block, model_language.Layer) else max(len(block.sequence), 1) for block in orso_blocks)
+    if total_layers < 2:
         raise ValueError('ORSO stack must contain at least 2 layers (superphase and subphase).')
 
-    logger.debug(f'Resolved layers: {orso_layers}')
-
-    # Convert ORSO layers to EasyReflectometry layers
-    erl_layers = []
-    for layer in orso_layers:
-        erl_layer = _convert_orso_layer_to_erl(layer)
-        erl_layers.append(erl_layer)
+    logger.debug(f'Resolved blocks: {orso_blocks}')
 
     # Create Superphase from first layer (thickness=0, roughness=0, both fixed)
-    superphase_layer = erl_layers[0]
+    superphase_layer = _convert_orso_layer_to_erl(orso_blocks[0])
     superphase_layer.thickness.value = 0.0
     superphase_layer.roughness.value = 0.0
     superphase_layer.thickness.fixed = True
@@ -135,51 +500,139 @@ def load_orso_model(orso_data) -> Sample:
     superphase = Multilayer(superphase_layer, name='Superphase')
 
     # Create Subphase from last layer (thickness=0 fixed, roughness enabled)
-    subphase_layer = erl_layers[-1]
+    subphase_layer = _convert_orso_layer_to_erl(orso_blocks[-1])
     subphase_layer.thickness.value = 0.0
     subphase_layer.thickness.fixed = True
     subphase_layer.roughness.fixed = False
     subphase = Multilayer(subphase_layer, name='Subphase')
 
+    # Middle blocks: consecutive plain layers group into one Multilayer;
+    # sub-stacks become their own (repeating) assemblies.
+    middle_assemblies = []
+    pending_layers = []
+
+    def flush_pending():
+        if pending_layers:
+            middle_assemblies.append(Multilayer(list(pending_layers), name='Loaded layer'))
+            pending_layers.clear()
+
+    for block in orso_blocks[1:-1]:
+        if isinstance(block, model_language.Layer):
+            pending_layers.append(_convert_orso_layer_to_erl(block))
+        else:
+            flush_pending()
+            middle_assemblies.append(_convert_orso_substack_to_erl(block))
+    flush_pending()
+
+    # Keep the historic single-group name; disambiguate multiple plain groups.
+    plain_groups = [a for a in middle_assemblies if a.name == 'Loaded layer']
+    if len(plain_groups) > 1:
+        for k, assembly in enumerate(plain_groups):
+            assembly.name = f'Loaded layer {k}'
+
     # Create Sample from the file
     sample_info = orso_data[0].info.data_source.sample
     sample_name = sample_info.name if sample_info.name else 'ORSO Sample'
 
-    # Build Sample based on number of layers
-    if len(erl_layers) == 2:
-        # Only superphase and subphase, no middle layers
-        sample = Sample(superphase, subphase, name=sample_name)
-    else:
-        # Create middle layer assembly from layers between first and last
-        middle_layers = erl_layers[1:-1]
-        loaded_layer = Multilayer(middle_layers, name='Loaded layer')
-        sample = Sample(superphase, loaded_layer, subphase, name=sample_name)
-
+    sample = Sample(superphase, *middle_assemblies, subphase, name=sample_name)
     return sample
 
 
+def _generate_layer_material(orso_layer) -> None:
+    """Resolve a model-language layer's material in place (formula -> SLD).
+
+    This is the per-layer half of orsopy's ``resolve_to_layers()``; a failure
+    is warned about (the material then imports with its declared values, or an
+    SLD of 0) instead of silently re-resolving the stack a different way.
+    """
+    try:
+        if orso_layer.material is None:
+            orso_layer.generate_material()
+        orso_layer.material.generate_density()
+    except Exception as e:
+        warnings.warn(
+            f'Could not resolve material for ORSO layer {getattr(orso_layer, "original_name", None) or orso_layer!r}: {e}',
+            UserWarning,
+            stacklevel=3,
+        )
+
+
+def _convert_orso_substack_to_erl(block) -> Multilayer:
+    """Convert an ORSO SubStack block into a (repeating) multilayer assembly."""
+    repetitions = int(getattr(block, 'repetitions', 1) or 1)
+    # SubStack.resolve_to_layers() returns layers * repetitions; resolve one
+    # period by temporarily neutralizing the repetition count.
+    original_repetitions = block.repetitions
+    try:
+        block.repetitions = 1
+        period = block.resolve_to_layers()
+    except Exception as e:
+        warnings.warn(
+            f'Could not fully resolve ORSO sub-stack ({e}); using its raw layer sequence.',
+            UserWarning,
+            stacklevel=3,
+        )
+        period = [item for item in block.sequence if isinstance(item, model_language.Layer)]
+    finally:
+        block.repetitions = original_repetitions
+    erl_layers = [_convert_orso_layer_to_erl(orso_layer) for orso_layer in period]
+    name = getattr(block, 'original_name', None) or 'Multilayer'
+    if repetitions > 1:
+        return RepeatingMultilayer(erl_layers, repetitions=repetitions, name=name)
+    return Multilayer(erl_layers, name=name)
+
+
 def _convert_orso_layer_to_erl(layer):
-    r"""Helper function to convert an ORSO layer to an EasyReflectometry laye."""
+    r"""Helper function to convert an ORSO layer to an EasyReflectometry layer."""
     material = layer.material
     # Prefer original_name for the material name, fall back to the formula; a
     # material defined only by its SLD has neither, so never leave it None.
-    m_name = layer.original_name if layer.original_name is not None else material.formula
+    formula = getattr(material, 'formula', None)
+    m_name = layer.original_name if layer.original_name is not None else formula
     if m_name is None:
         m_name = 'material'
 
-    # Get SLD values (use formula for density calculation if available)
-    formula_for_calc = material.formula if material.formula is not None else m_name
-    m_sld, m_isld = _get_sld_values(material, formula_for_calc)
+    erl_material = _convert_orso_material_to_erl(material, m_name)
 
-    # Create and return ERL layer
+    # Create and return ERL layer; lengths honour the declared unit (nm default).
     return Layer(
-        material=Material(sld=m_sld, isld=m_isld, name=m_name),
-        thickness=layer.thickness.magnitude if layer.thickness is not None else 0.0,
-        roughness=layer.roughness.magnitude if layer.roughness is not None else 0.0,
+        material=erl_material,
+        thickness=_length_to_angstrom(layer.thickness),
+        roughness=_length_to_angstrom(layer.roughness),
         name=layer.original_name if layer.original_name is not None else m_name,
     )
 
 
+def _convert_orso_material_to_erl(material, material_name):
+    """Convert an ORSO material to an ERL material.
+
+    Density-defined materials (formula + mass density, no SLD) stay
+    density-defined as :class:`MaterialDensity`, so formula and density remain
+    editable/recoverable instead of being flattened to a numeric SLD.
+    """
+    formula = getattr(material, 'formula', None)
+    mass_density = getattr(material, 'mass_density', None)
+    if material.sld is None and mass_density is not None and (formula or material_name):
+        magnitude = getattr(mass_density, 'magnitude', mass_density)
+        unit = getattr(mass_density, 'unit', None)
+        try:
+            factor = _MASS_DENSITY_UNIT_TO_G_CM3[unit]
+        except KeyError:
+            warnings.warn(
+                f'Unknown ORSO mass density unit {unit!r}; value used as g/cm^3.',
+                UserWarning,
+                stacklevel=4,
+            )
+            factor = 1.0
+        return MaterialDensity(
+            chemical_structure=formula if formula is not None else material_name,
+            density=float(magnitude) * factor,
+            name=material_name if material_name is not None else formula,
+        )
+    m_sld, m_isld = _get_sld_values(material, material_name)
+    return Material(sld=m_sld, isld=m_isld, name=material_name)
+
+
 def _get_sld_values(material, material_name):
     """Extract SLD values from material, calculating from density if needed
 
@@ -190,7 +643,7 @@ def _get_sld_values(material, material_name):
     if material.sld is None and material.mass_density is not None:
         # Calculate SLD from mass density
         # MaterialDensity already returns values in 10^-6 A^-2 units
-        m_density = material.mass_density.magnitude
+        m_density = getattr(material.mass_density, 'magnitude', material.mass_density)
         density = MaterialDensity(chemical_structure=material_name, density=m_density)
         m_sld = density.sld.value
         m_isld = density.isld.value
@@ -199,16 +652,26 @@ def _get_sld_values(material, material_name):
         m_sld = 0.0
         m_isld = 0.0
     else:
-        # ORSO stores SLD in absolute units (A^-2)
-        # Convert to internal representation (10^-6 A^-2) by multiplying by 1e6
-        if isinstance(material.sld, ComplexValue):
-            raw_sld = material.sld.real
-            m_sld = raw_sld * 1e6
-            m_isld = (material.sld.imag or 0.0) * 1e6
+        # ORSO stores SLD in absolute units (A^-2, or 1/nm^2 when declared).
+        # Convert to internal representation (10^-6 A^-2).
+        sld = material.sld
+        unit = getattr(sld, 'unit', None)
+        try:
+            unit_factor = _SLD_UNIT_TO_INV_ANGSTROM_SQUARED[unit]
+        except KeyError:
+            warnings.warn(
+                f'Unknown ORSO SLD unit {unit!r}; value used as 1/angstrom^2.',
+                UserWarning,
+                stacklevel=4,
+            )
+            unit_factor = 1.0
+        if isinstance(sld, ComplexValue):
+            raw_sld = sld.real
+            m_sld = raw_sld * unit_factor * 1e6
+            m_isld = (sld.imag if sld.imag is not None else 0.0) * unit_factor * 1e6
         else:
-            # A plain number, or an orsopy ``Value`` (unwrap its magnitude).
-            raw_sld = getattr(material.sld, 'magnitude', material.sld)
-            m_sld = raw_sld * 1e6
+            raw_sld = getattr(sld, 'magnitude', sld)
+            m_sld = raw_sld * unit_factor * 1e6
             m_isld = 0.0
         if raw_sld != 0.0 and abs(raw_sld) > 1e-2:
             warnings.warn(
@@ -222,45 +685,207 @@ def _get_sld_values(material, material_name):
     return m_sld, m_isld
 
 
-def load_orso_data(orso_data) -> DataSet1D:
-    """Convert parsed ORSO dataset objects into a scipp DataGroup.
+# ---------------------------------------------------------------------------
+# Sample -> model language, and .ort/.orb export
+# ---------------------------------------------------------------------------
+
+
+def _sanitize_model_key(name: str) -> str:
+    """Make a layer name safe for use in the model-language stack string."""
+    cleaned = ''.join('_' if character in '|()' else character for character in str(name)).strip()
+    return cleaned or 'layer'
+
+
+def _unique_model_key(existing: dict, name: str) -> str:
+    """A stack key not yet present in *existing* (appends a counter if needed)."""
+    key = _sanitize_model_key(name)
+    if key not in existing:
+        return key
+    counter = 2
+    while f'{key}_{counter}' in existing:
+        counter += 1
+    return f'{key}_{counter}'
+
+
+def _material_to_orso(material) -> model_language.Material:
+    """Convert an ERL material to an ORSO model-language material."""
+    if isinstance(material, MaterialDensity):
+        return model_language.Material(
+            formula=material.chemical_structure,
+            mass_density=Value(float(material.density.value), 'g/cm^3'),
+        )
+    return model_language.Material(
+        sld=ComplexValue(
+            real=float(material.sld.value) * 1e-6,
+            imag=float(material.isld.value) * 1e-6,
+            unit='1/angstrom^2',
+        )
+    )
+
+
+def sample_to_orso_model(sample: Sample) -> model_language.SampleModel:
+    """Convert an ERL Sample into an ORSO simple-model ``SampleModel`` (slab model).
+
+    Lengths are written in angstrom (declared via ``globals.length_unit``), SLDs
+    in 1/angstrom^2.  :class:`RepeatingMultilayer` assemblies are written with
+    the inline repetition syntax ``N ( layer1 | layer2 )``.
 
     Parameters
     ----------
-    orso_data : list
-        Parsed ORSO dataset list (as returned by ``orso.load_orso``).
+    sample : Sample
+        The sample to convert.
 
     Returns
     -------
-    sc.DataGroup
-        A scipp DataGroup with data, coords, and attrs.
+    model_language.SampleModel
+        The ORSO model-language representation.
     """
-    data = {}
-    coords = {}
-    attrs = {}
-    for i, o in enumerate(orso_data):
-        name = i
-        if o.info.data_set is not None:
-            name = o.info.data_set
-        coords[f'Qz_{name}'] = sc.array(
-            dims=[f'{o.info.columns[0].name}_{name}'],
-            values=o.data[:, 0],
-            variances=np.square(o.data[:, 3]),
-            unit=sc.Unit(o.info.columns[0].unit),
+    layers = {}
+    stack_parts = []
+
+    def add_layer(erl_layer) -> str:
+        key = _unique_model_key(layers, erl_layer.name or erl_layer.material.name)
+        layers[key] = model_language.Layer(
+            thickness=Value(float(erl_layer.thickness.value), 'angstrom'),
+            roughness=Value(float(erl_layer.roughness.value), 'angstrom'),
+            material=_material_to_orso(erl_layer.material),
         )
+        return key
+
+    for assembly in sample:
+        keys = [add_layer(erl_layer) for erl_layer in assembly.layers]
+        if isinstance(assembly, RepeatingMultilayer) and int(assembly.repetitions.value) > 1:
+            stack_parts.append(f'{int(assembly.repetitions.value)} ( ' + ' | '.join(keys) + ' )')
+        else:
+            stack_parts.extend(keys)
+
+    return model_language.SampleModel(
+        stack=' | '.join(stack_parts),
+        layers=layers,
+        globals=model_language.ModelParameters(length_unit='angstrom', sld_unit='1/angstrom^2'),
+    )
+
+
+_STANDARD_COLUMNS = [
+    Column(name='Qz', unit='1/angstrom', physical_quantity='wavevector transfer'),
+    Column(name='R', physical_quantity='reflectivity'),
+    ErrorColumn(error_of='R', error_type='uncertainty', value_is='sigma'),
+    ErrorColumn(error_of='Qz', error_type='resolution', value_is='sigma'),
+]
+
+
+def _orso_header_for_dataset(dataset: DataSet1D, model=None) -> Orso:
+    """Build an Orso header for one dataset, reusing a preserved header when present.
+
+    A header captured at load time (``dataset.orso_header``) keeps
+    ``data_source``/``reduction`` provenance; otherwise a minimal header is
+    synthesized.  Columns are always (re)set to the standard four with the
+    sigma convention, matching the exported data.
+    """
+    import copy as _copy
+
+    header_dict = getattr(dataset, 'orso_header', None)
+    if header_dict:
         try:
-            data[f'R_{name}'] = sc.array(
-                dims=[f'{o.info.columns[0].name}_{name}'],
-                values=o.data[:, 1],
-                variances=np.square(o.data[:, 2]),
-                unit=sc.Unit(o.info.columns[1].unit),
-            )
-        except TypeError:
-            data[f'R_{name}'] = sc.array(
-                dims=[f'{o.info.columns[0].name}_{name}'],
-                values=o.data[:, 1],
-                variances=np.square(o.data[:, 2]),
-            )
-        attrs[f'R_{name}'] = {'orso_header': sc.scalar(Header.asdict(o.info))}
-    data_group = sc.DataGroup(data=data, coords=coords, attrs=attrs)
-    return data_group
+            info = Orso.from_dict(_copy.deepcopy(dict(header_dict)))
+        except Exception:
+            info = Orso.empty()
+    else:
+        info = Orso.empty()
+    info.columns = _copy.deepcopy(_STANDARD_COLUMNS)
+    if model is not None and getattr(model, 'sample', None) is not None:
+        info.data_source.sample.model = sample_to_orso_model(model.sample)
+        if info.data_source.sample.name is None:
+            info.data_source.sample.name = model.sample.name
+    return info
+
+
+def _dataset_to_data_array(dataset: DataSet1D) -> np.ndarray:
+    """Assemble the Qz/R/sR/sQz array for export.
+
+    ``DataSet1D`` stores **variances** in ``ye``/``xe``; ORSO columns are sigma,
+    so the square root is taken here.  Absent errors are written as nan
+    (the spec's "unknown" marker), never as zeros.
+    """
+    x = np.asarray(dataset.x, dtype=float)
+    y = np.asarray(dataset.y, dtype=float)
+
+    def sigma_column(variances) -> np.ndarray:
+        if variances is None:
+            return np.full_like(x, np.nan)
+        variances = np.asarray(variances, dtype=float)
+        if variances.size != x.size or not np.any(np.nan_to_num(variances) > 0):
+            return np.full_like(x, np.nan)
+        return np.sqrt(variances)
+
+    return np.column_stack([x, y, sigma_column(dataset.ye), sigma_column(dataset.xe)])
+
+
+def orso_datasets_from_experiment(experiment, model=None) -> List[OrsoDataset]:
+    """Convert an experiment into a list of ``OrsoDataset`` objects.
+
+    A plain :class:`DataSet1D` becomes one dataset.  A ``PolarizedDataSet``
+    becomes one dataset per spin channel in a single multi-dataset file --
+    the format's intended packing for spin states -- with ``data_set`` labels
+    and per-dataset ``instrument_settings.polarization`` set to the channel.
+
+    Parameters
+    ----------
+    experiment :
+        A DataSet1D or PolarizedDataSet.
+    model : optional
+        The model whose sample is written into ``sample.model``. Defaults to
+        ``experiment.model``.
+
+    Returns
+    -------
+    List[OrsoDataset]
+        The datasets ready for ``save_orso``/``save_nexus``.
+    """
+    from orsopy.fileio.data_source import InstrumentSettings
+    from orsopy.fileio.data_source import Polarization
+
+    from easyreflectometry.data.polarized import PolarizedDataSet
+
+    if model is None:
+        model = getattr(experiment, 'model', None)
+
+    datasets = []
+    if isinstance(experiment, PolarizedDataSet):
+        for channel, channel_dataset in experiment.channels.items():
+            info = _orso_header_for_dataset(channel_dataset, model)
+            info.data_set = channel.value
+            measurement = info.data_source.measurement
+            if measurement.instrument_settings is None:
+                measurement.instrument_settings = InstrumentSettings(incident_angle=None, wavelength=None)
+            measurement.instrument_settings.polarization = Polarization(channel.value)
+            datasets.append(OrsoDataset(info=info, data=_dataset_to_data_array(channel_dataset)))
+    else:
+        info = _orso_header_for_dataset(experiment, model)
+        if info.data_set is None:
+            info.data_set = 0
+        datasets.append(OrsoDataset(info=info, data=_dataset_to_data_array(experiment)))
+    return datasets
+
+
+def save_orso_experiment(experiment, fname: str, model=None) -> None:
+    """Write an experiment to an ORSO file (`.ort` text, or `.orb` binary).
+
+    Parameters
+    ----------
+    experiment :
+        A DataSet1D or PolarizedDataSet (the latter is written as one file with
+        one ``data_set`` block per spin channel).
+    fname : str
+        Destination path; a ``.orb`` extension selects the binary (NeXus/HDF5)
+        representation, anything else the text one.
+    model : optional
+        The model whose sample is exported as the ORSO ``sample.model``.
+        Defaults to the experiment's model.
+    """
+    fname = str(fname)
+    datasets = orso_datasets_from_experiment(experiment, model=model)
+    if fname.lower().endswith('.orb'):
+        orso.save_nexus(datasets, fname)
+    else:
+        orso.save_orso(datasets, fname)
diff --git a/src/easyreflectometry/project.py b/src/easyreflectometry/project.py
index e5250a91..21ef26ee 100644
--- a/src/easyreflectometry/project.py
+++ b/src/easyreflectometry/project.py
@@ -30,8 +30,11 @@
 from easyreflectometry.data import DataSet1D
 from easyreflectometry.data import PolarizedDataSet
 from easyreflectometry.data import detect_polarization_channel
+from easyreflectometry.data import detect_polarization_channels_per_dataset
 from easyreflectometry.data import load_as_dataset
+from easyreflectometry.data.measurement import dataset_from_datagroup
 from easyreflectometry.data.measurement import extract_orso_title
+from easyreflectometry.data.measurement import load as load_measurement_file
 from easyreflectometry.data.measurement import load_data_from_orso_file
 from easyreflectometry.fitting import MultiFitter
 from easyreflectometry.inequality_constraints import InequalityEvaluation
@@ -44,6 +47,8 @@
 from easyreflectometry.model import ModelCollection
 from easyreflectometry.model import PercentageFwhm
 from easyreflectometry.model import Pointwise
+from easyreflectometry.orso_utils import is_orso_file
+from easyreflectometry.orso_utils import save_orso_experiment
 from easyreflectometry.sample import Layer
 from easyreflectometry.sample import Material
 from easyreflectometry.sample import MaterialCollection
@@ -774,22 +779,36 @@ def get_index_d2o(self) -> int:
         return [material.name for material in self._materials].index('D2O')
 
     def load_orso_file(self, path: Union[Path, str]) -> None:
-        """Load an ORSO file and optionally create a model and a data from it."""
-        from easyreflectometry.orso_utils import LoadOrso
-
-        model, data = LoadOrso(path)
-        if model is not None:
-            if isinstance(model, Sample):
-                model = Model(sample=model, name=model.name)
-            self.models = ModelCollection([model])
+        """Load an ORSO file, creating a model from its ``sample.model`` (when present) and an experiment.
+
+        .. deprecated::
+            Use :meth:`load_new_experiment` (data) together with
+            :meth:`set_sample_from_orso` /
+            :func:`easyreflectometry.orso_utils.load_orso_model` (model)
+            instead. This wrapper now routes through the same
+            ``DataSet1D`` + title + resolution path as every other importer.
+
+        Parameters
+        ----------
+        path : Union[Path, str]
+            Path to the ORSO file.
+        """
+        warnings.warn(
+            'Project.load_orso_file is deprecated; use load_new_experiment (data) and '
+            'set_sample_from_orso/load_orso_model (model) instead.',
+            DeprecationWarning,
+            stacklevel=2,
+        )
+        from easyreflectometry.orso_utils import _load_orso_any
+        from easyreflectometry.orso_utils import load_orso_model
+
+        orso_data = _load_orso_any(str(path))
+        sample = load_orso_model(orso_data)
+        if sample is not None:
+            self.models = ModelCollection([Model(sample=sample, name=sample.name)])
         else:
             self.default_model()
-        if data is not None:
-            self._experiments[0] = data
-            self._experiments[0].name = 'Experiment from ORSO'
-            self._experiments[0].model = self.models[0]
-            self._with_experiments = True
-        pass
+        self.load_experiment_for_model_at_index(path, 0)
 
     def set_sample_from_orso(self, sample: Sample) -> None:
         """Replace the current project model collection with a single model built from an ORSO-parsed sample.
@@ -938,7 +957,9 @@ def _apply_resolution_function(
         model : Model
             The model whose resolution function is set.
         """
-        if experiment.xe is not None and np.any(experiment.xe):
+        # nan-robust gate: nan is truthy for np.any, but a nan-carrying xe must
+        # not build a Pointwise (np.interp would propagate the nan everywhere).
+        if experiment.xe is not None and np.any(np.nan_to_num(experiment.xe) > 0):
             model.resolution_function = Pointwise(q_data_points=[experiment.x, experiment.y, experiment.xe])
         else:
             model.resolution_function = PercentageFwhm(5.0)
@@ -949,16 +970,27 @@ def _auto_set_background(experiment: DataSet1D) -> None:
         if experiment.model is not None and len(experiment.y) > 0:
             experiment.model.background = max(np.min(experiment.y), 1e-10)
 
-    def load_new_experiment(self, path: Union[Path, str]) -> None:
-        """Load new experiment."""
-        new_experiment = load_as_dataset(str(path))
+    def load_new_experiment(self, path: Union[Path, str], data_group=None) -> None:
+        """Load new experiment.
+
+        Parameters
+        ----------
+        path : Union[Path, str]
+            Path to the experiment data file.
+        data_group :
+            Pre-loaded scipp DataGroup for *path* (avoids re-parsing the
+            file). By default, None.
+        """
+        if data_group is None:
+            data_group = load_measurement_file(str(path))
+        new_experiment = load_as_dataset(str(path), data_group=data_group)
         new_index = len(self._experiments)
 
         model_index = 0
         if new_index < len(self.models):
             model_index = new_index
 
-        self._apply_experiment_metadata(path, new_experiment, f'Experiment {new_index}')
+        self._apply_experiment_metadata(path, new_experiment, f'Experiment {new_index}', data_group=data_group)
         new_experiment.model = self.models[model_index]
         self._auto_set_background(new_experiment)
         self._experiments[new_index] = new_experiment
@@ -976,12 +1008,16 @@ def count_datasets_in_file(self, path: Union[Path, str]) -> int:
         Returns
         -------
         int
-            Number of datasets found; 1 if the file cannot be introspected.
+            Number of datasets found; 1 if a non-ORSO file cannot be
+            introspected. A corrupt ORSO file (banner present but unparsable)
+            raises instead of being silently miscounted.
         """
         try:
-            data_group = load_data_from_orso_file(str(path))
+            data_group = load_measurement_file(str(path))
             return len(data_group['data'])
         except Exception:
+            if is_orso_file(str(path)):
+                raise
             return 1
 
     def load_all_experiments_from_file(self, path: Union[Path, str]) -> int:
@@ -989,8 +1025,9 @@ def load_all_experiments_from_file(self, path: Union[Path, str]) -> int:
 
         For a multi-dataset ORSO file (e.g. a multi-angle measurement), each dataset is
         registered as an independent experiment.  All experiments share the model that is
-        currently selected.  Falls back to :meth:`load_new_experiment` for single-dataset
-        files or on any loading error.
+        currently selected.  Single-dataset files go through
+        :meth:`load_new_experiment`.  A corrupt ORSO file raises (no silent
+        plain-text fallback).
 
         Parameters
         ----------
@@ -1002,32 +1039,19 @@ def load_all_experiments_from_file(self, path: Union[Path, str]) -> int:
         int
             Number of experiments that were added.
         """
-        try:
-            data_group = load_data_from_orso_file(str(path))
-        except Exception:
-            self.load_new_experiment(path)
-            return 1
+        data_group = load_measurement_file(str(path))
 
         data_keys = sorted(data_group['data'].keys())
         if len(data_keys) <= 1:
-            self.load_new_experiment(path)
+            self.load_new_experiment(path, data_group=data_group)
             return 1
 
         model_index = self._current_model_index
         for data_key in data_keys:
-            coord_key = data_key.replace('R_', 'Qz_')
             new_index = len(self._experiments)
 
-            d = data_group['data'][data_key]
-            c = data_group['coords'][coord_key]
-
-            new_experiment = DataSet1D(
-                name=f'Experiment {new_index}',
-                x=c.values,
-                y=d.values,
-                ye=d.variances,
-                xe=c.variances if c.variances is not None else None,
-            )
+            new_experiment = dataset_from_datagroup(data_group, data_key=data_key)
+            new_experiment.name = f'Experiment {new_index}'
             self._apply_experiment_metadata(
                 path,
                 new_experiment,
@@ -1091,19 +1115,119 @@ def load_polarized_experiment(
         """
         paths = {PolarizationChannel(channel): path for channel, path in paths.items()}
         channels = {}
+        title_data_group = None
         for channel, path in paths.items():
+            # Parse each file once: dataset count, data, and title all come
+            # from the same DataGroup (previously up to 3-4 parses per file).
+            data_group = load_measurement_file(str(path))
             # One file per channel means one dataset per file; a multi-dataset ORSO
             # file has no defined channel-to-dataset assignment here.
-            if self.count_datasets_in_file(path) > 1:
+            if len(data_group['data']) > 1:
                 raise ValueError(
-                    f"File '{path}' contains multiple datasets; polarized loading requires one dataset "
-                    f'per channel file. Export the {channel.value} channel to its own file.'
+                    f"File '{path}' contains multiple datasets; use load_polarized_experiment_from_file "
+                    f'for a single multi-dataset file, or export the {channel.value} channel to its own file.'
                 )
-            dataset = load_as_dataset(str(path))
+            dataset = load_as_dataset(str(path), data_group=data_group)
             # Keep the source file visible per channel (file → channel provenance).
             dataset.name = f'{channel.value}: {Path(path).name}'
             channels[channel] = dataset
+            if title_data_group is None:
+                title_data_group = data_group
+
+        first_path = next(iter(paths.values()))
+        return self._register_polarized_experiment(channels, model_index, str(first_path), title_data_group)
+
+    def load_polarized_experiment_from_file(
+        self,
+        path: Union[Path, str],
+        model_index: Optional[int] = None,
+    ) -> int:
+        """Load a polarized experiment from a single multi-dataset ORSO file.
+
+        Multi-dataset packing is the ORSO format's intended way to store spin
+        states: each ``data_set:`` block is classified by its **own** header
+        (``instrument_settings.polarization``), so per-dataset overrides are
+        honoured. Only the fully-analysed cross-sections ``pp/pm/mp/mm`` are
+        mapped; ``po/mo/op/om/unpolarized/vector`` are not coerced.
+
+        Parameters
+        ----------
+        path : Union[Path, str]
+            Path to the multi-dataset ORSO file.
+        model_index : Optional[int], optional
+            Index of the model the experiment belongs to. By default, the
+            current model.
+
+        Returns
+        -------
+        int
+            Index of the newly loaded experiment.
+
+        Raises
+        ------
+        ValueError :
+            If any dataset lacks a mappable polarization header, or two
+            datasets declare the same channel. Such files cannot be classified
+            unambiguously; use :meth:`load_all_experiments_from_file` (e.g. for
+            multi-angle files) or per-channel files instead.
+        """
+        from easyreflectometry.orso_utils import _load_orso_any
+        from easyreflectometry.orso_utils import _orso_dataset_key
+        from easyreflectometry.orso_utils import load_orso_data
+
+        orso_data = _load_orso_any(str(path))
+        classified = detect_polarization_channels_per_dataset(orso_data)
+
+        channels = {}
+        for i, (declared, channel) in enumerate(classified):
+            label = _orso_dataset_key(orso_data[i], i)
+            if channel is None:
+                reason = 'declares no mappable spin channel' if declared else 'declares no polarization'
+                raise ValueError(
+                    f"Dataset '{label}' in '{path}' {reason} "
+                    f'(only pp/pm/mp/mm are mapped; po/mo/op/om/unpolarized are not coerced). '
+                    f'Use load_all_experiments_from_file for non-polarized multi-dataset files, '
+                    f'or assign channels explicitly via load_polarized_experiment.'
+                )
+            if channel in channels:
+                raise ValueError(
+                    f"Duplicate spin channel '{channel.value}' in '{path}': more than one dataset "
+                    f'declares it. Assign channels explicitly via load_polarized_experiment.'
+                )
+            channels[channel] = label
+
+        data_group = load_orso_data(orso_data)
+        channel_datasets = {}
+        for channel, label in channels.items():
+            dataset = dataset_from_datagroup(data_group, data_key=f'R_{label}')
+            dataset.name = f'{channel.value}: {label}'
+            channel_datasets[channel] = dataset
+
+        return self._register_polarized_experiment(channel_datasets, model_index, str(path), data_group)
+
+    def _register_polarized_experiment(self, channels, model_index, title_path, title_data_group) -> int:
+        """Shared tail of the polarized loaders: build, name, and register the experiment.
 
+        Background and resolution follow the first (in canonical order)
+        channel; per-channel resolution functions are not supported (one per
+        experiment).
+
+        Parameters
+        ----------
+        channels :
+            Channel → DataSet1D mapping.
+        model_index :
+            Index of the model, or None for the current one.
+        title_path :
+            Path used for the fallback experiment title lookup.
+        title_data_group :
+            Pre-loaded DataGroup for the title lookup (avoids re-parsing).
+
+        Returns
+        -------
+        int
+            Index of the newly loaded experiment.
+        """
         new_index = len(self._experiments)
         if model_index is None:
             model_index = self._current_model_index
@@ -1115,13 +1239,11 @@ def load_polarized_experiment(
             model=model,
         )
         # Name from the ORSO title of the first file, when available.
-        first_channel = experiment.available_channels[0]
-        first_path = paths[first_channel]
-        self._apply_experiment_metadata(first_path, experiment, f'Polarized experiment {new_index}')
+        self._apply_experiment_metadata(
+            title_path, experiment, f'Polarized experiment {new_index}', data_group=title_data_group
+        )
 
-        # Background and resolution follow the first (in canonical order) channel;
-        # per-channel resolution functions are not supported (one per experiment).
-        first_dataset = experiment[first_channel]
+        first_dataset = experiment[experiment.available_channels[0]]
         self._auto_set_background(first_dataset)
         self._apply_resolution_function(first_dataset, model)
 
@@ -1129,11 +1251,32 @@ def load_polarized_experiment(
         self._with_experiments = True
         return new_index
 
+    def save_experiment_as_orso(self, path: Union[Path, str], index: Optional[int] = None) -> None:
+        """Write an experiment to an ORSO file (`.ort` text, or `.orb` binary).
+
+        The experiment's model (when set) is exported as the ORSO
+        ``sample.model`` (slab representation); data columns are written as
+        sigma. A polarized experiment becomes one file with one ``data_set:``
+        block per spin channel.
+
+        Parameters
+        ----------
+        path : Union[Path, str]
+            Destination path; a ``.orb`` extension selects the binary format.
+        index : Optional[int], optional
+            Index of the experiment to save. By default, the current one.
+        """
+        if index is None:
+            index = self._current_experiment_index
+        experiment = self._experiments[index]
+        save_orso_experiment(experiment, str(path), model=experiment.model)
+
     def load_experiment_for_model_at_index(self, path: Union[Path, str], index: Optional[int] = 0) -> None:
         """Load experiment for model at index."""
-        experiment = load_as_dataset(str(path))
+        data_group = load_measurement_file(str(path))
+        experiment = load_as_dataset(str(path), data_group=data_group)
 
-        self._apply_experiment_metadata(path, experiment, f'Experiment {index}')
+        self._apply_experiment_metadata(path, experiment, f'Experiment {index}', data_group=data_group)
         experiment.model = self.models[index]
         self._auto_set_background(experiment)
         self._experiments[index] = experiment
diff --git a/tests/_static/example.orb b/tests/_static/example.orb
new file mode 100644
index 00000000..a4900027
Binary files /dev/null and b/tests/_static/example.orb differ
diff --git a/tests/_static/fwhm_sqz.ort b/tests/_static/fwhm_sqz.ort
new file mode 100644
index 00000000..e847b96f
--- /dev/null
+++ b/tests/_static/fwhm_sqz.ort
@@ -0,0 +1,39 @@
+# # ORSO reflectivity data file | 1.2 standard | YAML encoding | https://www.reflectometry.org/
+# data_source:
+#   owner:
+#     name: Generated fixture
+#     affiliation: EasyScience
+#   experiment:
+#     title: FWHM sQz fixture
+#     instrument: fixture
+#     start_date: null
+#     probe: neutron
+#   sample:
+#     name: FWHM sQz fixture
+#   measurement:
+#     instrument_settings:
+#       incident_angle: {magnitude: 0.5, unit: deg}
+#       wavelength: {magnitude: 4.0, unit: angstrom}
+#       polarization: null
+#     data_files: []
+# reduction:
+#   software: {name: make_fixtures}
+# data_set: 0
+# columns:
+# - {name: Qz, unit: 1/angstrom, physical_quantity: wavevector transfer}
+# - {name: R, physical_quantity: reflectivity}
+# - {error_of: R, error_type: uncertainty, value_is: sigma}
+# - {error_of: Qz, error_type: resolution, value_is: FWHM}
+# # Qz (1/angstrom)    R                      sR                     sQz                   
+1.0000000000000000e-02 6.7032104603563936e-01 3.3516052301781972e-02 4.7096400900618988e-04
+3.6363636363636362e-02 2.3350747909091338e-01 1.1675373954545669e-02 1.7125963963861449e-03
+6.2727272727272729e-02 8.1343153050479095e-02 4.0671576525239551e-03 2.9542287837661003e-03
+8.9090909090909082e-02 2.8336598603717900e-02 1.4168299301858950e-03 4.1958611711460552e-03
+1.1545454545454545e-01 9.8717265313318422e-03 4.9358632656659209e-04 5.4374935585260100e-03
+1.4181818181818182e-01 3.4394748181587544e-03 1.7197374090793772e-04 6.6791259459059657e-03
+1.6818181818181818e-01 1.1987952218190572e-03 5.9939761090952858e-05 7.9207583332859214e-03
+1.9454545454545455e-01 4.1825284298601604e-04 2.0912642149300804e-05 9.1623907206658754e-03
+2.2090909090909092e-01 1.4635033351987538e-04 7.3175166759937690e-06 1.0404023108045831e-02
+2.4727272727272728e-01 5.1632895160517905e-05 2.5816447580258954e-06 1.1645655495425785e-02
+2.7363636363636362e-01 1.8638006121158546e-05 9.3190030605792728e-07 1.2887287882805741e-02
+2.9999999999999999e-01 7.1442123533282095e-06 3.5721061766641050e-07 1.4128920270185696e-02
diff --git a/tests/_static/nan_sqz.ort b/tests/_static/nan_sqz.ort
new file mode 100644
index 00000000..c10b221d
--- /dev/null
+++ b/tests/_static/nan_sqz.ort
@@ -0,0 +1,39 @@
+# # ORSO reflectivity data file | 1.2 standard | YAML encoding | https://www.reflectometry.org/
+# data_source:
+#   owner:
+#     name: Generated fixture
+#     affiliation: EasyScience
+#   experiment:
+#     title: All-nan sQz fixture
+#     instrument: fixture
+#     start_date: null
+#     probe: neutron
+#   sample:
+#     name: All-nan sQz fixture
+#   measurement:
+#     instrument_settings:
+#       incident_angle: {magnitude: 0.5, unit: deg}
+#       wavelength: {magnitude: 4.0, unit: angstrom}
+#       polarization: null
+#     data_files: []
+# reduction:
+#   software: {name: make_fixtures}
+# data_set: 0
+# columns:
+# - {name: Qz, unit: 1/angstrom, physical_quantity: wavevector transfer}
+# - {name: R, physical_quantity: reflectivity}
+# - {error_of: R, error_type: uncertainty, value_is: sigma}
+# - {error_of: Qz, error_type: resolution, value_is: sigma}
+# # Qz (1/angstrom)    R                      sR                     sQz                   
+1.0000000000000000e-02 6.7032104603563936e-01 3.3516052301781972e-02 nan                   
+3.6363636363636362e-02 2.3350747909091338e-01 1.1675373954545669e-02 nan                   
+6.2727272727272729e-02 8.1343153050479095e-02 4.0671576525239551e-03 nan                   
+8.9090909090909082e-02 2.8336598603717900e-02 1.4168299301858950e-03 nan                   
+1.1545454545454545e-01 9.8717265313318422e-03 4.9358632656659209e-04 nan                   
+1.4181818181818182e-01 3.4394748181587544e-03 1.7197374090793772e-04 nan                   
+1.6818181818181818e-01 1.1987952218190572e-03 5.9939761090952858e-05 nan                   
+1.9454545454545455e-01 4.1825284298601604e-04 2.0912642149300804e-05 nan                   
+2.2090909090909092e-01 1.4635033351987538e-04 7.3175166759937690e-06 nan                   
+2.4727272727272728e-01 5.1632895160517905e-05 2.5816447580258954e-06 nan                   
+2.7363636363636362e-01 1.8638006121158546e-05 9.3190030605792728e-07 nan                   
+2.9999999999999999e-01 7.1442123533282095e-06 3.5721061766641050e-07 nan                   
diff --git a/tests/_static/nm_units.ort b/tests/_static/nm_units.ort
new file mode 100644
index 00000000..95a6f212
--- /dev/null
+++ b/tests/_static/nm_units.ort
@@ -0,0 +1,67 @@
+# # ORSO reflectivity data file | 1.2 standard | YAML encoding | https://www.reflectometry.org/
+# data_source:
+#   owner:
+#     name: Generated fixture
+#     affiliation: EasyScience
+#   experiment:
+#     title: nm units fixture
+#     instrument: fixture
+#     start_date: null
+#     probe: neutron
+#   sample:
+#     name: nm units fixture
+#     model:
+#       stack: air | 3 ( A | B ) | film | Si
+#       layers:
+#         air:
+#           thickness: 0.0
+#           roughness: 0.0
+#           material:
+#             sld: {real: 0.0, imag: 0.0}
+#         A:
+#           thickness: 2.0
+#           roughness: 0.3
+#           material:
+#             sld: {real: 4.0e-06, imag: 0.0}
+#         B:
+#           thickness: 1.0
+#           roughness: 0.3
+#           material:
+#             sld: {real: 2.0e-06, imag: 0.0}
+#         film:
+#           thickness: 10.0
+#           roughness: 0.5
+#           material:
+#             sld: {real: 3.47e-06, imag: 0.0}
+#         Si:
+#           thickness: 0.0
+#           roughness: 0.3
+#           material:
+#             sld: {real: 2.07e-06, imag: 0.0}
+#   measurement:
+#     instrument_settings:
+#       incident_angle: {magnitude: 0.5, unit: deg}
+#       wavelength: {magnitude: 4.0, unit: angstrom}
+#       polarization: null
+#     data_files: []
+# reduction:
+#   software: {name: make_fixtures}
+# data_set: 0
+# columns:
+# - {name: Qz, unit: 1/nm, physical_quantity: wavevector transfer}
+# - {name: R, physical_quantity: reflectivity}
+# - {error_of: R, error_type: uncertainty, value_is: sigma}
+# - {error_of: Qz, error_type: resolution, value_is: sigma}
+# # Qz (1/nm)          R                      sR                     sQz                   
+1.0000000000000001e-01 6.7032104603563936e-01 3.3516052301781972e-02 2.0000000000000000e-03
+3.6363636363636365e-01 2.3350747909091338e-01 1.1675373954545669e-02 7.2727272727272727e-03
+6.2727272727272732e-01 8.1343153050479095e-02 4.0671576525239551e-03 1.2545454545454545e-02
+8.9090909090909087e-01 2.8336598603717900e-02 1.4168299301858950e-03 1.7818181818181816e-02
+1.1545454545454545e+00 9.8717265313318422e-03 4.9358632656659209e-04 2.3090909090909092e-02
+1.4181818181818182e+00 3.4394748181587544e-03 1.7197374090793772e-04 2.8363636363636365e-02
+1.6818181818181819e+00 1.1987952218190572e-03 5.9939761090952858e-05 3.3636363636363638e-02
+1.9454545454545455e+00 4.1825284298601604e-04 2.0912642149300804e-05 3.8909090909090914e-02
+2.2090909090909090e+00 1.4635033351987538e-04 7.3175166759937690e-06 4.4181818181818183e-02
+2.4727272727272727e+00 5.1632895160517905e-05 2.5816447580258954e-06 4.9454545454545452e-02
+2.7363636363636363e+00 1.8638006121158546e-05 9.3190030605792728e-07 5.4727272727272722e-02
+3.0000000000000000e+00 7.1442123533282095e-06 3.5721061766641050e-07 5.9999999999999998e-02
diff --git a/tests/_static/partial_nan_sqz.ort b/tests/_static/partial_nan_sqz.ort
new file mode 100644
index 00000000..895ace5c
--- /dev/null
+++ b/tests/_static/partial_nan_sqz.ort
@@ -0,0 +1,39 @@
+# # ORSO reflectivity data file | 1.2 standard | YAML encoding | https://www.reflectometry.org/
+# data_source:
+#   owner:
+#     name: Generated fixture
+#     affiliation: EasyScience
+#   experiment:
+#     title: Partial-nan sQz fixture
+#     instrument: fixture
+#     start_date: null
+#     probe: neutron
+#   sample:
+#     name: Partial-nan sQz fixture
+#   measurement:
+#     instrument_settings:
+#       incident_angle: {magnitude: 0.5, unit: deg}
+#       wavelength: {magnitude: 4.0, unit: angstrom}
+#       polarization: null
+#     data_files: []
+# reduction:
+#   software: {name: make_fixtures}
+# data_set: 0
+# columns:
+# - {name: Qz, unit: 1/angstrom, physical_quantity: wavevector transfer}
+# - {name: R, physical_quantity: reflectivity}
+# - {error_of: R, error_type: uncertainty, value_is: sigma}
+# - {error_of: Qz, error_type: resolution, value_is: sigma}
+# # Qz (1/angstrom)    R                      sR                     sQz                   
+1.0000000000000000e-02 6.7032104603563936e-01 3.3516052301781972e-02 nan                   
+3.6363636363636362e-02 2.3350747909091338e-01 1.1675373954545669e-02 nan                   
+6.2727272727272729e-02 8.1343153050479095e-02 4.0671576525239551e-03 1.2545454545454546e-03
+8.9090909090909082e-02 2.8336598603717900e-02 1.4168299301858950e-03 1.7818181818181817e-03
+1.1545454545454545e-01 9.8717265313318422e-03 4.9358632656659209e-04 2.3090909090909091e-03
+1.4181818181818182e-01 3.4394748181587544e-03 1.7197374090793772e-04 2.8363636363636364e-03
+1.6818181818181818e-01 1.1987952218190572e-03 5.9939761090952858e-05 3.3636363636363638e-03
+1.9454545454545455e-01 4.1825284298601604e-04 2.0912642149300804e-05 3.8909090909090911e-03
+2.2090909090909092e-01 1.4635033351987538e-04 7.3175166759937690e-06 4.4181818181818185e-03
+2.4727272727272728e-01 5.1632895160517905e-05 2.5816447580258954e-06 4.9454545454545454e-03
+2.7363636363636362e-01 1.8638006121158546e-05 9.3190030605792728e-07 5.4727272727272723e-03
+2.9999999999999999e-01 7.1442123533282095e-06 3.5721061766641050e-07 nan                   
diff --git a/tests/_static/polarized_2ch.ort b/tests/_static/polarized_2ch.ort
new file mode 100644
index 00000000..fd91ffc8
--- /dev/null
+++ b/tests/_static/polarized_2ch.ort
@@ -0,0 +1,57 @@
+# # ORSO reflectivity data file | 1.2 standard | YAML encoding | https://www.reflectometry.org/
+# data_source:
+#   owner:
+#     name: Generated fixture
+#     affiliation: EasyScience
+#   experiment:
+#     title: Polarized fixture
+#     instrument: fixture
+#     start_date: null
+#     probe: neutron
+#   sample:
+#     name: Polarized fixture
+#   measurement:
+#     instrument_settings:
+#       incident_angle: {magnitude: 0.5, unit: deg}
+#       wavelength: {magnitude: 4.0, unit: angstrom}
+#       polarization: pp
+#     data_files: []
+# reduction:
+#   software: {name: make_fixtures}
+# data_set: pp
+# columns:
+# - {name: Qz, unit: 1/angstrom, physical_quantity: wavevector transfer}
+# - {name: R, physical_quantity: reflectivity}
+# - {error_of: R, error_type: uncertainty, value_is: sigma}
+# - {error_of: Qz, error_type: resolution, value_is: sigma}
+# # Qz (1/angstrom)    R                      sR                     sQz                   
+1.0000000000000000e-02 7.3735315063920337e-01 3.3516052301781972e-02 2.0000000000000001e-04
+3.6363636363636362e-02 2.5685822700000471e-01 1.1675373954545669e-02 7.2727272727272723e-04
+6.2727272727272729e-02 8.9477468355527015e-02 4.0671576525239551e-03 1.2545454545454546e-03
+8.9090909090909082e-02 3.1170258464089693e-02 1.4168299301858950e-03 1.7818181818181817e-03
+1.1545454545454545e-01 1.0858899184465027e-02 4.9358632656659209e-04 2.3090909090909091e-03
+1.4181818181818182e-01 3.7834222999746300e-03 1.7197374090793772e-04 2.8363636363636364e-03
+1.6818181818181818e-01 1.3186747440009630e-03 5.9939761090952858e-05 3.3636363636363638e-03
+1.9454545454545455e-01 4.6007812728461768e-04 2.0912642149300804e-05 3.8909090909090911e-03
+2.2090909090909092e-01 1.6098536687186293e-04 7.3175166759937690e-06 4.4181818181818185e-03
+2.4727272727272728e-01 5.6796184676569701e-05 2.5816447580258954e-06 4.9454545454545454e-03
+2.7363636363636362e-01 2.0501806733274400e-05 9.3190030605792728e-07 5.4727272727272723e-03
+2.9999999999999999e-01 7.8586335886610303e-06 3.5721061766641050e-07 6.0000000000000001e-03
+# data_set: mm
+# data_source:
+#   measurement:
+#     instrument_settings:
+#       polarization: mm
+# # Qz (1/angstrom)    R                      sR                     sQz                   
+1.0000000000000000e-02 6.0328894143207545e-01 3.3516052301781972e-02 2.0000000000000001e-04
+3.6363636363636362e-02 2.1015673118182204e-01 1.1675373954545669e-02 7.2727272727272723e-04
+6.2727272727272729e-02 7.3208837745431188e-02 4.0671576525239551e-03 1.2545454545454546e-03
+8.9090909090909082e-02 2.5502938743346110e-02 1.4168299301858950e-03 1.7818181818181817e-03
+1.1545454545454545e-01 8.8845538781986578e-03 4.9358632656659209e-04 2.3090909090909091e-03
+1.4181818181818182e-01 3.0955273363428792e-03 1.7197374090793772e-04 2.8363636363636364e-03
+1.6818181818181818e-01 1.0789156996371516e-03 5.9939761090952858e-05 3.3636363636363638e-03
+1.9454545454545455e-01 3.7642755868741445e-04 2.0912642149300804e-05 3.8909090909090911e-03
+2.2090909090909092e-01 1.3171530016788786e-04 7.3175166759937690e-06 4.4181818181818185e-03
+2.4727272727272728e-01 4.6469605644466116e-05 2.5816447580258954e-06 4.9454545454545454e-03
+2.7363636363636362e-01 1.6774205509042691e-05 9.3190030605792728e-07 5.4727272727272723e-03
+2.9999999999999999e-01 6.4297911179953888e-06 3.5721061766641050e-07 6.0000000000000001e-03
diff --git a/tests/test_orso_support.py b/tests/test_orso_support.py
new file mode 100644
index 00000000..aea94c01
--- /dev/null
+++ b/tests/test_orso_support.py
@@ -0,0 +1,339 @@
+# SPDX-FileCopyrightText: 2026 EasyScience contributors <https://github.com/easyscience>
+# SPDX-License-Identifier: BSD-3-Clause
+
+"""Tests for the ORSO support update (ORSO_UPDATE_TASK.md):
+
+reader hardening (banner discriminator, FWHM->sigma, nan policy, 1/nm Q),
+model-language preservation (units, repeats, density materials),
+single-file polarized import, the .ort/.orb exporter, and .orb reading.
+"""
+
+import os
+
+import numpy as np
+import pytest
+from easyscience import global_object
+
+import easyreflectometry
+from easyreflectometry.data import DataSet1D
+from easyreflectometry.data import PolarizedDataSet
+from easyreflectometry.data.measurement import load
+from easyreflectometry.data.measurement import load_as_dataset
+from easyreflectometry.model import PercentageFwhm
+from easyreflectometry.model import Pointwise
+from easyreflectometry.orso_utils import SIGMA_TO_FWHM
+from easyreflectometry.orso_utils import _load_orso_any
+from easyreflectometry.orso_utils import is_orso_file
+from easyreflectometry.orso_utils import load_orso_model
+from easyreflectometry.orso_utils import sample_to_orso_model
+from easyreflectometry.orso_utils import save_orso_experiment
+from easyreflectometry.project import Project
+from easyreflectometry.sample import Layer
+from easyreflectometry.sample import Material
+from easyreflectometry.sample import Multilayer
+from easyreflectometry.sample import RepeatingMultilayer
+from easyreflectometry.sample import Sample
+from easyreflectometry.sample.elements.materials.material_density import MaterialDensity
+
+PATH_STATIC = os.path.join(os.path.dirname(easyreflectometry.__file__), '..', '..', 'tests', '_static')
+
+# The Qz/sQz grids the generated fixtures were built from (see task 7a).
+FIXTURE_QZ = np.linspace(0.01, 0.3, 12)
+FIXTURE_SQZ = FIXTURE_QZ * 0.02
+
+
+@pytest.fixture(autouse=True)
+def clear_global_map():
+    global_object.map._clear()
+    yield
+    global_object.map._clear()
+
+
+@pytest.fixture
+def project() -> Project:
+    return Project()
+
+
+class TestBannerDiscriminator:
+    def test_ort_file_is_orso(self):
+        assert is_orso_file(os.path.join(PATH_STATIC, 'Ni_example.ort')) is True
+
+    def test_orb_file_is_orso(self):
+        assert is_orso_file(os.path.join(PATH_STATIC, 'example.orb')) is True
+
+    def test_txt_file_is_not_orso(self):
+        assert is_orso_file(os.path.join(PATH_STATIC, 'test_example1.txt')) is False
+
+    def test_bannered_but_corrupt_file_raises(self, tmp_path):
+        # A file carrying the ORSO banner that fails to parse must raise, not
+        # silently fall back to plain-text loading (which drops the header).
+        bad = tmp_path / 'bad.ort'
+        bad.write_text(
+            '# # ORSO reflectivity data file | 1.1 standard | YAML encoding | https://www.reflectometry.org/\n'
+            '# data_source: {[[not yaml\n'
+            '1.0 2.0 3.0 4.0\n'
+        )
+        with pytest.raises(ValueError, match='Error loading ORSO file'):
+            load(str(bad))
+
+    def test_non_bannered_file_falls_back_to_txt(self, tmp_path):
+        plain = tmp_path / 'plain.dat'
+        plain.write_text('0.01 1.0 0.1 0.001\n0.02 0.5 0.05 0.002\n')
+        data_group = load(str(plain))
+        assert 'R_plain' in data_group['data']
+
+    def test_count_datasets_raises_for_corrupt_orso(self, project, tmp_path):
+        bad = tmp_path / 'bad.ort'
+        bad.write_text('# # ORSO reflectivity data file | 1.1 standard | YAML encoding\n# data_source: {[[\n1 2 3 4\n')
+        with pytest.raises(ValueError):
+            project.count_datasets_in_file(str(bad))
+
+    def test_count_datasets_multi(self, project):
+        assert project.count_datasets_in_file(os.path.join(PATH_STATIC, 'polarized_2ch.ort')) == 2
+
+
+class TestNanPolicy:
+    def test_all_nan_sqz_leaves_xe_empty(self):
+        data_group = load(os.path.join(PATH_STATIC, 'nan_sqz.ort'))
+        coords = data_group['coords'][list(data_group['coords'])[0]]
+        assert coords.variances is None
+
+    def test_all_nan_sqz_falls_back_to_percentage_fwhm(self, project):
+        project.default_model()
+        project.load_new_experiment(os.path.join(PATH_STATIC, 'nan_sqz.ort'))
+        assert isinstance(project.models[0].resolution_function, PercentageFwhm)
+
+    def test_partial_nan_sqz_is_interpolated(self):
+        with pytest.warns(UserWarning, match='sQz values are nan'):
+            data_group = load(os.path.join(PATH_STATIC, 'partial_nan_sqz.ort'))
+        coords = data_group['coords'][list(data_group['coords'])[0]]
+        assert coords.variances is not None
+        assert np.all(np.isfinite(coords.variances))
+
+    def test_partial_nan_sqz_builds_finite_pointwise(self, project):
+        project.default_model()
+        with pytest.warns(UserWarning, match='sQz values are nan'):
+            project.load_new_experiment(os.path.join(PATH_STATIC, 'partial_nan_sqz.ort'))
+        resolution_function = project.models[0].resolution_function
+        assert isinstance(resolution_function, Pointwise)
+        assert np.all(np.isfinite(resolution_function.smearing(FIXTURE_QZ)))
+
+
+class TestValueIsFwhm:
+    def test_fwhm_sqz_is_converted_to_sigma(self):
+        data_group = load(os.path.join(PATH_STATIC, 'fwhm_sqz.ort'))
+        coords = data_group['coords'][list(data_group['coords'])[0]]
+        # The fixture stores sigma * SIGMA_TO_FWHM declared as FWHM; loading
+        # must convert back so stored variances are sigma squared.
+        np.testing.assert_allclose(np.sqrt(coords.variances), FIXTURE_SQZ)
+
+    def test_sigma_file_is_not_scaled(self):
+        data_group = load(os.path.join(PATH_STATIC, 'partial_nan_sqz.ort'))
+        coords = data_group['coords'][list(data_group['coords'])[0]]
+        valid = np.sqrt(coords.variances)[2:-1]
+        np.testing.assert_allclose(valid, FIXTURE_SQZ[2:-1])
+
+
+class TestQUnitConversion:
+    def test_1_per_nm_qz_is_converted_to_1_per_angstrom(self):
+        data_group = load(os.path.join(PATH_STATIC, 'nm_units.ort'))
+        coords = data_group['coords'][list(data_group['coords'])[0]]
+        np.testing.assert_allclose(coords.values, FIXTURE_QZ)
+        np.testing.assert_allclose(np.sqrt(coords.variances), FIXTURE_SQZ)
+        assert str(coords.unit) == '1/Å'
+
+
+class TestModelLanguagePreservation:
+    def test_nm_default_length_unit_is_honoured(self):
+        # nm_units.ort declares thicknesses as bare magnitudes with no
+        # length_unit -> the model-language default (nm) applies.
+        sample = load_orso_model(_load_orso_any(os.path.join(PATH_STATIC, 'nm_units.ort')))
+        film = [layer for assembly in sample for layer in assembly.layers if layer.name == 'film'][0]
+        assert film.thickness.value == pytest.approx(100.0)  # 10 nm
+        assert film.roughness.value == pytest.approx(5.0)  # 0.5 nm
+
+    def test_default_roughness_from_globals(self):
+        sample = load_orso_model(_load_orso_any(os.path.join(PATH_STATIC, 'nm_units.ort')))
+        subphase_layer = sample[-1].layers[0]
+        # Si declares roughness 0.3 (nm default) -> 3 A
+        assert subphase_layer.roughness.value == pytest.approx(3.0)
+
+    def test_repeated_substack_becomes_repeating_multilayer(self):
+        sample = load_orso_model(_load_orso_any(os.path.join(PATH_STATIC, 'nm_units.ort')))
+        repeating = [assembly for assembly in sample if isinstance(assembly, RepeatingMultilayer)]
+        assert len(repeating) == 1
+        assert repeating[0].repetitions.value == 3
+        assert [layer.name for layer in repeating[0].layers] == ['A', 'B']
+        assert repeating[0].layers[0].thickness.value == pytest.approx(20.0)  # 2 nm
+
+    def test_density_material_stays_density_defined(self):
+        sample = load_orso_model(_load_orso_any(os.path.join(PATH_STATIC, 'Ni_example.ort')))
+        m1 = sample[1].layers[0]
+        assert isinstance(m1.material, MaterialDensity)
+        assert m1.material.chemical_structure == 'Ni'
+        assert m1.material.density.value == pytest.approx(8.9)
+        # and the derived SLD is still sensible (Ni ~ 9.4e-6 A^-2)
+        assert m1.material.sld.value == pytest.approx(9.4, abs=0.1)
+
+
+class TestPolarizedSingleFile:
+    def test_load_polarized_experiment_from_file(self, project):
+        project.default_model()
+        index = project.load_polarized_experiment_from_file(os.path.join(PATH_STATIC, 'polarized_2ch.ort'))
+        experiment = project.experiments[index]
+        assert isinstance(experiment, PolarizedDataSet)
+        assert [channel.value for channel in experiment.available_channels] == ['pp', 'mm']
+        assert experiment.name == 'Polarized fixture'
+        assert experiment.model is project.models[0]
+        np.testing.assert_allclose(experiment['pp'].x, FIXTURE_QZ)
+
+    def test_unclassifiable_dataset_raises(self, project):
+        project.default_model()
+        # nan_sqz.ort carries no spin-channel polarization -> not coerced.
+        with pytest.raises(ValueError, match='declares no (mappable spin channel|polarization)'):
+            project.load_polarized_experiment_from_file(os.path.join(PATH_STATIC, 'nan_sqz.ort'))
+
+    def test_multidataset_file_rejected_by_per_channel_loader(self, project):
+        project.default_model()
+        with pytest.raises(ValueError, match='load_polarized_experiment_from_file'):
+            project.load_polarized_experiment({'pp': os.path.join(PATH_STATIC, 'polarized_2ch.ort')})
+
+    def test_duplicate_channel_raises(self, project, tmp_path):
+        # Build a file where two datasets declare the same channel.
+        datasets = _load_orso_any(os.path.join(PATH_STATIC, 'polarized_2ch.ort'))
+        from orsopy.fileio.data_source import Polarization
+        from orsopy.fileio.orso import save_orso
+
+        for orso_dataset in datasets:
+            orso_dataset.info.data_source.measurement.instrument_settings.polarization = Polarization('pp')
+        duplicate_file = tmp_path / 'dup.ort'
+        save_orso(datasets, str(duplicate_file))
+
+        project.default_model()
+        with pytest.raises(ValueError, match="Duplicate spin channel 'pp'"):
+            project.load_polarized_experiment_from_file(str(duplicate_file))
+
+
+class TestExporter:
+    def test_ort_data_roundtrip_writes_sigma(self, project, tmp_path):
+        project.default_model()
+        project.load_new_experiment(os.path.join(PATH_STATIC, 'Ni_example.ort'))
+        experiment = project.experiments[0]
+        out = tmp_path / 'out.ort'
+        project.save_experiment_as_orso(str(out), 0)
+
+        back = _load_orso_any(str(out))
+        assert len(back) == 1
+        data = back[0].data
+        np.testing.assert_allclose(data[:, 0], experiment.x)
+        np.testing.assert_allclose(data[:, 1], experiment.y)
+        # DataSet1D stores variances; the file must carry sigma.
+        np.testing.assert_allclose(data[:, 2], np.sqrt(experiment.ye))
+        np.testing.assert_allclose(data[:, 3], np.sqrt(experiment.xe))
+        # sigma convention declared in the columns
+        assert back[0].info.columns[2].value_is == 'sigma'
+        assert back[0].info.columns[3].value_is == 'sigma'
+
+    def test_export_reuses_preserved_header(self, project, tmp_path):
+        project.default_model()
+        project.load_new_experiment(os.path.join(PATH_STATIC, 'Ni_example.ort'))
+        out = tmp_path / 'out.ort'
+        project.save_experiment_as_orso(str(out), 0)
+        info = _load_orso_any(str(out))[0].info
+        # data_source/reduction provenance from the original file, not synthesized
+        assert info.data_source.experiment.title == 'Metal films'
+        assert info.data_source.owner.name == 'Joe Bloggs'
+
+    def test_export_writes_model_language(self, project, tmp_path):
+        project.default_model()
+        project.load_new_experiment(os.path.join(PATH_STATIC, 'Ni_example.ort'))
+        out = tmp_path / 'out.ort'
+        project.save_experiment_as_orso(str(out), 0)
+        model = _load_orso_any(str(out))[0].info.data_source.sample.model
+        assert model is not None
+        assert model.globals.length_unit == 'angstrom'
+        # the exported model resolves back into layers
+        assert len(model.resolve_to_layers()) >= 2
+
+    def test_absent_errors_written_as_nan(self, tmp_path):
+        dataset = DataSet1D(x=np.array([0.01, 0.02]), y=np.array([1.0, 0.5]))
+        out = tmp_path / 'nan.ort'
+        save_orso_experiment(dataset, str(out))
+        data = _load_orso_any(str(out))[0].data
+        assert np.all(np.isnan(data[:, 2]))
+        assert np.all(np.isnan(data[:, 3]))
+
+    def test_polarized_export_single_multidataset_file(self, project, tmp_path):
+        project.default_model()
+        index = project.load_polarized_experiment_from_file(os.path.join(PATH_STATIC, 'polarized_2ch.ort'))
+        out = tmp_path / 'pol.ort'
+        project.save_experiment_as_orso(str(out), index)
+
+        back = _load_orso_any(str(out))
+        assert [d.info.data_set for d in back] == ['pp', 'mm']
+        polarizations = [str(d.info.data_source.measurement.instrument_settings.polarization.value) for d in back]
+        assert polarizations == ['pp', 'mm']
+
+    def test_repeating_multilayer_roundtrips_via_model_language(self):
+        air = Layer(material=Material(sld=0.0, isld=0.0, name='air'), thickness=0, roughness=0, name='air')
+        layer_a = Layer(material=Material(sld=4.0, isld=0.0, name='A'), thickness=20, roughness=3, name='A')
+        layer_b = Layer(material=Material(sld=2.0, isld=0.0, name='B'), thickness=10, roughness=3, name='B')
+        si = Layer(material=Material(sld=2.07, isld=0.0, name='Si'), thickness=0, roughness=3, name='Si')
+        sample = Sample(
+            Multilayer(air, name='Superphase'),
+            RepeatingMultilayer([layer_a, layer_b], repetitions=5, name='rep'),
+            Multilayer(si, name='Subphase'),
+            name='test',
+        )
+        orso_model = sample_to_orso_model(sample)
+        assert '5 ( A | B )' in orso_model.stack
+        assert len(orso_model.resolve_to_layers()) == 12  # 1 + 5*2 + 1
+
+    def test_model_as_orso_returns_model_language_dict(self, project):
+        project.default_model()
+        orso_dict = project.models[0].as_orso()
+        assert 'stack' in orso_dict
+        assert 'layers' in orso_dict
+        assert orso_dict['globals']['length_unit'] == 'angstrom'
+
+
+class TestOrbSupport:
+    def test_load_orb_file(self):
+        data_group = load(os.path.join(PATH_STATIC, 'example.orb'))
+        assert 'R_0' in data_group['data']
+        coords = data_group['coords'][list(data_group['coords'])[0]]
+        np.testing.assert_allclose(coords.values, FIXTURE_QZ)
+
+    def test_orb_write_and_read_roundtrip(self, project, tmp_path):
+        project.default_model()
+        project.load_new_experiment(os.path.join(PATH_STATIC, 'Ni_example.ort'))
+        experiment = project.experiments[0]
+        out = tmp_path / 'out.orb'
+        project.save_experiment_as_orso(str(out), 0)
+        back = _load_orso_any(str(out))
+        np.testing.assert_allclose(back[0].data[:, 0], experiment.x)
+        np.testing.assert_allclose(back[0].data[:, 2], np.sqrt(experiment.ye))
+
+    def test_load_orb_as_dataset(self):
+        dataset = load_as_dataset(os.path.join(PATH_STATIC, 'example.orb'))
+        assert isinstance(dataset, DataSet1D)
+        np.testing.assert_allclose(dataset.x, FIXTURE_QZ)
+        assert dataset.orso_header is not None
+
+
+class TestFwhmSigmaContract:
+    def test_pointwise_serialization_stays_variances(self):
+        # The saved-project contract: sQz_data_points round-trip as variances.
+        qz = np.array([0.01, 0.02])
+        reflectivity = np.array([1.0, 0.5])
+        variances = np.array([1e-8, 2e-8])
+        pointwise = Pointwise([qz, reflectivity, variances])
+        as_dict = pointwise.as_dict()
+        np.testing.assert_allclose(as_dict['sQz_data_points'], variances)
+        from easyreflectometry.model.resolution_functions import ResolutionFunction
+
+        restored = ResolutionFunction.from_dict(as_dict)
+        np.testing.assert_allclose(restored.smearing(qz), np.sqrt(variances))
+
+    def test_sigma_to_fwhm_constant(self):
+        assert SIGMA_TO_FWHM == pytest.approx(2.3548, abs=1e-4)
diff --git a/tests/unit/test_project_core.py b/tests/unit/test_project_core.py
index 5a227255..91383dbc 100644
--- a/tests/unit/test_project_core.py
+++ b/tests/unit/test_project_core.py
@@ -105,15 +105,31 @@ def test_percentage_fwhm_fallback_when_q_variances_absent(self, project: Project
 
 class TestLoadOrsoFile:
     def test_load_orso_file_creates_model_and_experiment(self, project: Project):
-        with pytest.warns(UserWarning):
+        # example.ort has no sample.model -> default model; deprecated wrapper
+        # now goes through the DataSet1D + title + resolution path.
+        with pytest.warns((UserWarning, DeprecationWarning)):
             project.load_orso_file(os.path.join(PATH_STATIC, 'example.ort'))
 
         assert len(project.models) == 1
         assert len(project.experiments) == 1
-        assert project.experiments[0].name == 'Experiment from ORSO'
+        assert isinstance(project.experiments[0], DataSet1D)
+        assert project.experiments[0].name == 'Example data file from refnx docs'
         assert project.experiments[0].model is project.models[0]
         assert project._with_experiments is True
 
+    def test_load_orso_file_with_model_builds_sample_from_file(self, project: Project):
+        # Ni_example.ort carries a sample.model -> the model comes from the file.
+        with pytest.warns(DeprecationWarning):
+            project.load_orso_file(os.path.join(PATH_STATIC, 'Ni_example.ort'))
+
+        assert len(project.models) == 1
+        assert project.models[0].sample.name == 'Ni on Si'
+        assert len(project.experiments) == 1
+        assert isinstance(project.experiments[0], DataSet1D)
+        assert project.experiments[0].model is project.models[0]
+        # The measured sQz feeds a Pointwise resolution function.
+        assert isinstance(project.models[0].resolution_function, Pointwise)
+
 
 class TestModelData:
     def test_model_data_for_model_at_index_returns_reflectivity(self, project: Project):