diff --git a/docs/user-guide/estia-data-reduction.md b/docs/user-guide/estia-data-reduction.md
new file mode 100644
index 00000000..35e429f0
--- /dev/null
+++ b/docs/user-guide/estia-data-reduction.md
@@ -0,0 +1,321 @@
+# Polarization data reduction for ESTIA
+
+Based on https://confluence.ess.eu/display/ESTIA/Polarised+Neutron+Reflectometry+%28PNR%29+-+Reduction+Notes
+
+## Model
+
+Intensity in the detector is related to the reflectivity of the sample by the model
+```math
+\begin{bmatrix}
+I^{+} \\
+I^{-}
+\end{bmatrix}
+\big(\lambda, j\big)
+=
+D(\lambda, j)
+\begin{bmatrix}
+1 - a^{\uparrow} & 1 - a^{\downarrow} \\
+a^{\uparrow} & a^{\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+(1 - f_2) & f_2 \\ f_2 & (1 - f_2)
+\end{bmatrix}
+\begin{bmatrix}
+R^{\uparrow\uparrow} & R^{\downarrow\uparrow} \\
+R^{\uparrow\downarrow} & R^{\downarrow\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+(1 - f_1) & f_1 \\ f_1 & (1 - f_1)
+\end{bmatrix}
+\begin{bmatrix}
+1 - p^{\uparrow} \\
+1 - p^{\downarrow}
+\end{bmatrix}
+```
+where
+
+* $I^+$ is the intensity of the neutron beam transmitted by the analyzer
+and $I^-$ is the intensity of the neutron beam reflected by the analyzer,
+* $R^\cdot$ are the reflectivities of the sample,
+  - $R^{\uparrow\uparrow}$ is the fraction of incoming neutrons with spin up that are reflected with spin up,
+  - $R^{\uparrow\downarrow}$ is the fraction of incoming neutrons with spin up that are reflected with spin down,
+  - etc..
+* $a^\uparrow$ is the analyzer reflectivity for spin up neutrons and $a^\downarrow$ is the analyzer reflectivity for spin down neutrons,
+* $p^\uparrow$ is the polarizer reflectivity for spin up neutrons and $p^\downarrow$ is the polarizer reflectivity for spin down neutrons,
+* $f_1$ is the probability of spin flip by the polarizer spin flipper, $f_2$ is the probability of spin flip by the analyzer spin flipper
+* $D$ represents the inhomogeneity from the beam- and detector efficiency (and all other polarization unrelated terms).
+
+## Reducing a measurement
+
+If the sample is measured at two different flipper settings $f_1=0, f_2=0$ and $f_1=1, f_2=0$, then we have four measurement in total:
+```math
+\begin{bmatrix}
+I^{0+} \\
+I^{0-}
+\end{bmatrix}
+\big(\lambda, j\big)
+=
+D(\lambda, j)
+\begin{bmatrix}
+1 - a^{\uparrow} & 1 - a^{\downarrow} \\
+a^{\uparrow} & a^{\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+R^{\uparrow\uparrow} & R^{\downarrow\uparrow} \\
+R^{\uparrow\downarrow} & R^{\downarrow\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+1 - p^{\uparrow} \\
+1 - p^{\downarrow}
+\end{bmatrix}
+```
+```math
+\begin{bmatrix}
+I^{1+} \\
+I^{1-}
+\end{bmatrix}
+\big(\lambda, j\big)
+=
+D(\lambda, j)
+\begin{bmatrix}
+1 - a^{\uparrow} & 1 - a^{\downarrow} \\
+a^{\uparrow} & a^{\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+R^{\uparrow\uparrow} & R^{\downarrow\uparrow} \\
+R^{\uparrow\downarrow} & R^{\downarrow\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+1 - p^{\downarrow} \\
+1 - p^{\uparrow}
+\end{bmatrix}.
+```
+
+To simplify the above, collect the terms in the matrix $\mathbf{a}$
+```math
+\begin{bmatrix}
+I^{0+} \\
+I^{0-} \\
+I^{1+} \\
+I^{1-}
+\end{bmatrix}
+\big(\lambda, j\big)
+=
+D(\lambda, j)
+\mathbf{a}(\lambda)
+\begin{bmatrix}
+R^{\uparrow\uparrow} \\
+R^{\uparrow\downarrow} \\
+R^{\downarrow\uparrow} \\
+R^{\downarrow\downarrow}
+\end{bmatrix}
+(Q(\lambda, j)).
+```
+
+To compute the reflectivities, integrate over a region of (almost) constant $Q$
+```math
+\int_{Q\in[q_{n}, q_{n+1}]}
+\mathbf{a}^{-1}(\lambda)
+\begin{bmatrix}
+I^{0+} \\
+I^{0-} \\
+I^{1+} \\
+I^{1-}
+\end{bmatrix}
+\big(\lambda, j\big)
+d\lambda dj
+\approx
+\int_{Q\in[q_{n}, q_{n+1}]}
+D(\lambda, j)
+d\lambda dj
+\begin{bmatrix}
+R^{\uparrow\uparrow} \\
+R^{\uparrow\downarrow} \\
+R^{\downarrow\uparrow} \\
+R^{\downarrow\downarrow}
+\end{bmatrix}
+(q_{n+\frac{1}{2}}).
+```
+The integral on the righ-hand-side can be evaluated using the reference measurement, call evaluated integral $\bar{D}(q_{n+{\frac{1}{2}}})$.
+$R$ was moved outside of the integral because if $Q$ is almost constant so is $R(Q)$.
+
+Finally we have
+```math
+\int_{Q\in[q_{n}, q_{n+1}]}
+\mathbf{a}^{-1}(\lambda)
+\bar{D}^{-1}(q_{n+{\frac{1}{2}}})
+\begin{bmatrix}
+I^{0+} \\
+I^{0-} \\
+I^{1+} \\
+I^{1-}
+\end{bmatrix}
+\big(\lambda, j\big)
+d\lambda dj
+\approx
+\begin{bmatrix}
+R^{\uparrow\uparrow} \\
+R^{\uparrow\downarrow} \\
+R^{\downarrow\uparrow} \\
+R^{\downarrow\downarrow}
+\end{bmatrix}
+(q_{n+\frac{1}{2}}).
+```
+
+### How to use the reference measurement to compute the integral over $D(\lambda, j)$?
+
+For a reference measurement using flipper setting $f_1=0, f_2=0$ we have
+```math
+\begin{bmatrix}
+I_{ref}^{+} \\
+I_{ref}^{-}
+\end{bmatrix}
+\big(\lambda, j\big)
+=
+D(\lambda, j)
+\begin{bmatrix}
+1 - a^{\uparrow} & 1 - a^{\downarrow} \\
+a^{\uparrow} & a^{\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+R_{ref}^{\uparrow\uparrow} & R_{ref}^{\downarrow\uparrow} \\
+R_{ref}^{\uparrow\downarrow} & R_{ref}^{\downarrow\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+1 - p^{\uparrow} \\
+1 - p^{\downarrow}
+\end{bmatrix}.
+```
+But in practice, the analyzer/polarizer will be efficient enough to make only one of $I_{ref}^\pm$ have enough intensity to be useful. For example:
+```math
+\frac{I_{ref}^{+}(\lambda, j)}{r^+(\lambda, j)}
+=
+D(\lambda, j)
+```
+where $r^+$ is a known term involving the reflectivity of the supermirror and the pol-/analyzer efficiencies.
+The expression for $D$ above can be used to evaluate integrals of $D$,
+but in this case only in the region of the detector where the transmitted beam hits, because we only got data in that region from our reference measurement.
+
+To measure $D$ for the entire detector we need to make several reference measurements with different flipper settings so that every part of the detector is illuminated in at least one measurement.
+It might be unecessary to use all 4 flipper settings, but to illustrate the idea imagine we make reference measurements using all 4 flipper settings:
+```math
+\begin{bmatrix}
+I_{ref}^{00+} \\
+I_{ref}^{00-}
+\end{bmatrix}
+\big(\lambda, j\big)
+=
+D(\lambda, j)
+\begin{bmatrix}
+1 - a^{\uparrow} & 1 - a^{\downarrow} \\
+a^{\uparrow} & a^{\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+R_{ref}^{\uparrow\uparrow} & R_{ref}^{\downarrow\uparrow} \\
+R_{ref}^{\uparrow\downarrow} & R_{ref}^{\downarrow\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+1 - p^{\uparrow} \\
+1 - p^{\downarrow}
+\end{bmatrix}
+```
+```math
+\begin{bmatrix}
+I_{ref}^{01+} \\
+I_{ref}^{01-}
+\end{bmatrix}
+\big(\lambda, j\big)
+=
+D(\lambda, j)
+\begin{bmatrix}
+1 - a^{\uparrow} & 1 - a^{\downarrow} \\
+a^{\uparrow} & a^{\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+0 & 1 \\ 1 & 0
+\end{bmatrix}
+\begin{bmatrix}
+R_{ref}^{\uparrow\uparrow} & R_{ref}^{\downarrow\uparrow} \\
+R_{ref}^{\uparrow\downarrow} & R_{ref}^{\downarrow\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+1 - p^{\uparrow} \\
+1 - p^{\downarrow}
+\end{bmatrix}
+```
+```math
+\begin{bmatrix}
+I_{ref}^{10+} \\
+I_{ref}^{10-}
+\end{bmatrix}
+\big(\lambda, j\big)
+=
+D(\lambda, j)
+\begin{bmatrix}
+1 - a^{\uparrow} & 1 - a^{\downarrow} \\
+a^{\uparrow} & a^{\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+R_{ref}^{\uparrow\uparrow} & R_{ref}^{\downarrow\uparrow} \\
+R_{ref}^{\uparrow\downarrow} & R_{ref}^{\downarrow\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+0 & 1 \\ 1 & 0
+\end{bmatrix}
+\begin{bmatrix}
+1 - p^{\uparrow} \\
+1 - p^{\downarrow}
+\end{bmatrix}
+```
+```math
+\begin{bmatrix}
+I_{ref}^{11+} \\
+I_{ref}^{11-}
+\end{bmatrix}
+\big(\lambda, j\big)
+=
+D(\lambda, j)
+\begin{bmatrix}
+1 - a^{\uparrow} & 1 - a^{\downarrow} \\
+a^{\uparrow} & a^{\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+0 & 1 \\ 1 & 0
+\end{bmatrix}
+\begin{bmatrix}
+R_{ref}^{\uparrow\uparrow} & R_{ref}^{\downarrow\uparrow} \\
+R_{ref}^{\uparrow\downarrow} & R_{ref}^{\downarrow\downarrow}
+\end{bmatrix}
+\begin{bmatrix}
+0 & 1 \\ 1 & 0
+\end{bmatrix}
+\begin{bmatrix}
+1 - p^{\uparrow} \\
+1 - p^{\downarrow}
+\end{bmatrix}.
+```
+
+Summing all 8 measurements gives us an expression for $D$ that ought to be valid for the entire detector:
+```math
+\frac{
+I_{ref}^{00+}(\lambda, j) +
+I_{ref}^{00-}(\lambda, j) +
+I_{ref}^{01+}(\lambda, j) +
+I_{ref}^{01-}(\lambda, j) +
+I_{ref}^{10+}(\lambda, j) +
+I_{ref}^{10-}(\lambda, j) +
+I_{ref}^{11+}(\lambda, j) +
+I_{ref}^{11-}(\lambda, j)
+}{
+r^{00+}(\lambda, j) +
+r^{00-}(\lambda, j) +
+r^{01+}(\lambda, j) +
+r^{01-}(\lambda, j) +
+r^{10+}(\lambda, j) +
+r^{10-}(\lambda, j) +
+r^{11+}(\lambda, j) +
+r^{11-}(\lambda, j)
+}
+=
+D(\lambda, j).
+```
diff --git a/src/ess/estia/__init__.py b/src/ess/estia/__init__.py
new file mode 100644
index 00000000..dc33be71
--- /dev/null
+++ b/src/ess/estia/__init__.py
@@ -0,0 +1,85 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
+import importlib.metadata
+
+import sciline
+import scipp as sc
+
+from ..reflectometry import providers as reflectometry_providers
+from ..reflectometry import supermirror
+from ..reflectometry.types import (
+    BeamSize,
+    DetectorSpatialResolution,
+    NeXusDetectorName,
+    RunType,
+    SamplePosition,
+    BeamDivergenceLimits,
+)
+from . import conversions, load, orso, resolution, utils, figures
+from .instrument_view import instrument_view
+from .types import (
+    AngularResolution,
+    SampleSizeResolution,
+    WavelengthResolution,
+)
+
+
+try:
+    __version__ = importlib.metadata.version(__package__ or __name__)
+except importlib.metadata.PackageNotFoundError:
+    __version__ = "0.0.0"
+
+
+providers = (
+    *reflectometry_providers,
+    *load.providers,
+    *conversions.providers,
+    *resolution.providers,
+    *utils.providers,
+    *figures.providers,
+    *orso.providers,
+)
+"""
+List of providers for setting up a Sciline pipeline.
+
+This provides a default Estia workflow including providers for loadings files.
+"""
+
+
+def default_parameters() -> dict:
+    return {
+        supermirror.MValue: sc.scalar(5, unit=sc.units.dimensionless),
+        supermirror.CriticalEdge: 0.022 * sc.Unit("1/angstrom"),
+        supermirror.Alpha: sc.scalar(0.25 / 0.088, unit=sc.units.angstrom),
+        BeamSize[RunType]: 2.0 * sc.units.mm,
+        DetectorSpatialResolution[RunType]: 0.0025 * sc.units.m,
+        SamplePosition[RunType]: sc.vector([0, 0, 0], unit="m"),
+        NeXusDetectorName[RunType]: "detector",
+        BeamDivergenceLimits: (
+            sc.scalar(-0.75, unit='deg'),
+            sc.scalar(0.75, unit='deg'),
+        ),
+    }
+
+
+def EstiaWorkflow() -> sciline.Pipeline:
+    """
+    Workflow with default parameters for the Estia instrument.
+    """
+    return sciline.Pipeline(providers=providers, params=default_parameters())
+
+
+__all__ = [
+    "supermirror",
+    "conversions",
+    "load",
+    "orso",
+    "resolution",
+    "instrument_view",
+    "providers",
+    "default_parameters",
+    "WavelengthResolution",
+    "AngularResolution",
+    "SampleSizeResolution",
+    "EstiaWorkflow",
+]
diff --git a/src/ess/estia/conversions.py b/src/ess/estia/conversions.py
new file mode 100644
index 00000000..f843bed6
--- /dev/null
+++ b/src/ess/estia/conversions.py
@@ -0,0 +1,97 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
+import scipp as sc
+
+from ..reflectometry.conversions import reflectometry_q
+from ..reflectometry.types import (
+    BeamDivergenceLimits,
+    BeamSize,
+    RawDetectorData,
+    ReducibleData,
+    RunType,
+    WavelengthBins,
+    YIndexLimits,
+    ZIndexLimits,
+)
+
+
+def theta(divergence_angle, sample_rotation, detector_rotation):
+    '''
+    Angle of reflection.
+
+    Computes the angle between the scattering direction of
+    the neutron and the sample surface.
+    '''
+    return divergence_angle + detector_rotation - sample_rotation
+
+
+def angle_of_divergence(
+    theta, sample_rotation, angle_to_center_of_beam, natural_incidence_angle
+):
+    """
+    Difference between the incident angle and the center of the incident beam.
+    Useful for filtering parts of the beam that have too high divergence.
+
+    This is always in the interval [-0.75 deg, 0.75 deg],
+    but the divergence of the incident beam can also be reduced.
+    """
+    return (
+        theta
+        - sample_rotation
+        - angle_to_center_of_beam
+        - natural_incidence_angle.to(unit='rad')
+    )
+
+
+def wavelength(
+    event_time_offset,
+    # Other inputs
+):
+    "Converts event_time_offset to wavelength"
+    # Use frame unwrapping from scippneutron
+    pass
+
+
+def _not_between(v, a, b):
+    return (v < a) | (v > b)
+
+
+def add_common_coords_and_masks(
+    da: RawDetectorData[RunType],
+    ylim: YIndexLimits,
+    zlims: ZIndexLimits,
+    bdlim: BeamDivergenceLimits,
+    wbins: WavelengthBins,
+    beam_size: BeamSize[RunType],
+) -> ReducibleData[RunType]:
+    "Adds coords and masks that are useful for both reference and sample measurements."
+    da = da.transform_coords(
+        ("wavelength", "theta", "angle_of_divergence", "Q"),
+        {
+            "divergence_angle": "pixel_divergence_angle",
+            "wavelength": wavelength,
+            "theta": theta,
+            "angle_of_divergence": angle_of_divergence,
+            "Q": reflectometry_q,
+        },
+        rename_dims=False,
+        keep_intermediate=False,
+    )
+    da.masks["stripe_range"] = _not_between(da.coords["stripe"], *ylim)
+    da.masks['z_range'] = _not_between(da.coords["z_index"], *zlims)
+    da.bins.masks["divergence_too_large"] = _not_between(
+        da.bins.coords["angle_of_divergence"],
+        bdlim[0].to(unit=da.bins.coords["angle_of_divergence"].bins.unit),
+        bdlim[1].to(unit=da.bins.coords["angle_of_divergence"].bins.unit),
+    )
+    da.bins.masks['wavelength'] = _not_between(
+        da.bins.coords['wavelength'],
+        wbins[0],
+        wbins[-1],
+    )
+    # Correct for illumination of virtual source
+    da /= sc.sin(da.bins.coords['theta'])
+    return da
+
+
+providers = (add_common_coords_and_masks,)
diff --git a/src/ess/estia/normalization.py b/src/ess/estia/normalization.py
new file mode 100644
index 00000000..5ecaf455
--- /dev/null
+++ b/src/ess/estia/normalization.py
@@ -0,0 +1,183 @@
+import scipp as sc
+
+from ..reflectometry.conversions import reflectometry_q
+from ..reflectometry.supermirror import (
+    Alpha,
+    CriticalEdge,
+    MValue,
+    supermirror_reflectivity,
+)
+from ..reflectometry.types import (
+    DetectorSpatialResolution,
+    QBins,
+    ReducedReference,
+    ReducibleData,
+    Reference,
+    ReferenceRun,
+    ReflectivityOverQ,
+    ReflectivityOverZW,
+    Sample,
+    SampleRun,
+    WavelengthBins,
+)
+from .conversions import theta
+from .resolution import (
+    angular_resolution,
+    q_resolution,
+    sample_size_resolution,
+    wavelength_resolution,
+)
+
+
+def mask_events_where_supermirror_does_not_cover(
+    sam: ReducibleData[SampleRun],
+    ref: ReducedReference,
+    critical_edge: CriticalEdge,
+    mvalue: MValue,
+    alpha: Alpha,
+) -> Sample:
+    """
+    Mask events in regions of the detector the reference does not cover.
+
+    Regions of the detector that the reference
+    measurement doesn't cover cannot be used to compute reflectivity.
+
+    Preferably the reference measurement should cover the entire
+    detector, but sometimes that is not possible, for example
+    if the supermirror :math:`M` value was too limited or because the reference
+    was measured at too high angle.
+
+    To figure out what events need to be masked,
+    compute the supermirror reflectivity as a function
+    of the :math:`Q` the event would have had if it had belonged to
+    the reference measurement.
+    """
+    R = supermirror_reflectivity(
+        reflectometry_q(
+            sam.bins.coords["wavelength"],
+            theta(
+                sam.coords["pixel_divergence_angle"],
+                ref.coords["sample_rotation"],
+                ref.coords["detector_rotation"],
+            ),
+        ),
+        c=critical_edge,
+        m=mvalue,
+        alpha=alpha,
+    )
+    sam.bins.masks["supermirror_does_not_cover"] = sc.isnan(R)
+    return sam
+
+
+def reduce_reference(
+    reference: ReducibleData[ReferenceRun],
+    wavelength_bins: WavelengthBins,
+    critical_edge: CriticalEdge,
+    mvalue: MValue,
+    alpha: Alpha,
+) -> ReducedReference:
+    """
+    Reduces the reference measurement to estimate the
+    intensity distribution in the detector for
+    an ideal sample with reflectivity :math:`R = 1`.
+    """
+    R = supermirror_reflectivity(
+        reference.bins.coords['Q'],
+        c=critical_edge,
+        m=mvalue,
+        alpha=alpha,
+    )
+    reference.bins.masks['invalid'] = sc.isnan(R)
+    reference /= R
+    return reference.bins.concat(('stripe',)).hist(wavelength=wavelength_bins)
+
+
+def reduce_sample_over_q(
+    sample: Sample,
+    reference: Reference,
+    qbins: QBins,
+) -> ReflectivityOverQ:
+    """
+    Computes reflectivity as ratio of
+    sample intensity and intensity from a sample
+    with ideal reflectivity.
+
+    Returns reflectivity as a function of :math:`Q`.
+    """
+    h = reference.flatten(to='Q').hist(Q=qbins)
+    R = sample.bins.concat().bin(Q=qbins) / h.data
+    R.coords['Q_resolution'] = sc.sqrt(
+        (
+            (reference * reference.coords['Q_resolution'] ** 2)
+            .flatten(to='Q')
+            .hist(Q=qbins)
+        )
+        / h
+    ).data
+    return R
+
+
+def reduce_sample_over_zw(
+    sample: Sample,
+    reference: Reference,
+    wbins: WavelengthBins,
+) -> ReflectivityOverZW:
+    """
+    Computes reflectivity as ratio of
+    sample intensity and intensity from a sample
+    with ideal reflectivity.
+
+    Returns reflectivity as a function of ``blade``, ``wire`` and :math:`\\wavelength`.
+    """
+    R = sample.bins.concat(('stripe',)).bin(wavelength=wbins) / reference.data
+    R.masks["too_few_events"] = reference.data < sc.scalar(1, unit="counts")
+    return R
+
+
+def evaluate_reference(
+    reference: ReducedReference,
+    sample: ReducibleData[SampleRun],
+    qbins: QBins,
+    detector_spatial_resolution: DetectorSpatialResolution[SampleRun],
+) -> Reference:
+    """
+    Adds a :math:`Q` and :math:`Q`-resolution coordinate to each bin of the ideal
+    intensity distribution. The coordinates are computed as if the data came from
+    the sample measurement, that is, they use the ``sample_rotation``
+    and ``detector_rotation`` parameters from the sample measurement.
+    """
+    ref = reference.copy()
+    ref.coords["sample_rotation"] = sample.coords["sample_rotation"]
+    ref.coords["detector_rotation"] = sample.coords["detector_rotation"]
+    ref.coords["sample_size"] = sample.coords["sample_size"]
+    ref.coords["detector_spatial_resolution"] = detector_spatial_resolution
+    ref.coords["wavelength"] = sc.midpoints(ref.coords["wavelength"])
+    ref = ref.transform_coords(
+        (
+            "Q",
+            "wavelength_resolution",
+            "sample_size_resolution",
+            "angular_resolution",
+            "Q_resolution",
+        ),
+        {
+            "divergence_angle": "pixel_divergence_angle",
+            "theta": theta,
+            "Q": reflectometry_q,
+            "wavelength_resolution": wavelength_resolution,
+            "sample_size_resolution": sample_size_resolution,
+            "angular_resolution": angular_resolution,
+            "Q_resolution": q_resolution,
+        },
+        rename_dims=False,
+    )
+    return sc.values(ref)
+
+
+providers = (
+    reduce_reference,
+    reduce_sample_over_q,
+    reduce_sample_over_zw,
+    evaluate_reference,
+    mask_events_where_supermirror_does_not_cover,
+)
diff --git a/src/ess/estia/resolution.py b/src/ess/estia/resolution.py
new file mode 100644
index 00000000..3940e00f
--- /dev/null
+++ b/src/ess/estia/resolution.py
@@ -0,0 +1,116 @@
+# SPDX-License-Identifier: BSD-3-Clause
+# Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
+import scipp as sc
+
+from ..reflectometry.tools import fwhm_to_std
+
+
+def wavelength_resolution(
+    # What parameters are needed?
+):
+    """
+    Find the wavelength resolution contribution of the ESTIA instrument.
+
+    Parameters
+    ----------
+
+    L1:
+        Distance from midpoint between choppers to sample.
+    L2:
+        Distance from sample to detector.
+
+    Returns
+    -------
+    :
+        The wavelength resolution variable, as standard deviation.
+    """
+    # Don't yet know how to do this
+    raise NotImplementedError()
+
+
+def sample_size_resolution(
+    L2,
+    sample_size,
+):
+    """
+    The resolution from the projected sample size, where it may be bigger
+    than the detector pixel resolution as described in Section 4.3.3 of the Amor
+    publication (doi: 10.1016/j.nima.2016.03.007).
+
+    Parameters
+    ----------
+    L2:
+        Distance from sample to detector.
+    sample_size:
+        Size of sample.
+
+    Returns
+    -------
+    :
+        Standard deviation of contribution from the sample size.
+    """
+    return fwhm_to_std(sample_size / L2.to(unit=sample_size.unit))
+
+
+def angular_resolution(
+    theta,
+    L2,
+    detector_spatial_resolution,
+):
+    """
+    Determine the angular resolution of the ESTIA instrument.
+
+    Parameters
+    ----------
+    theta:
+        Angle of reflection.
+    L2:
+        Distance between sample and detector.
+    detector_spatial_resolution:
+        FWHM of detector pixel resolution.
+
+    Returns
+    -------
+    :
+        Angular resolution standard deviation
+    """
+    return (
+        fwhm_to_std(
+            sc.atan(
+                detector_spatial_resolution
+                / L2.to(unit=detector_spatial_resolution.unit)
+            )
+        ).to(unit=theta.unit)
+        / theta
+    )
+
+
+def q_resolution(
+    Q,
+    angular_resolution,
+    wavelength_resolution,
+    sample_size_resolution,
+):
+    """
+    Compute resolution in Q.
+
+    Parameters
+    ----------
+    Q:
+        Momentum transfer.
+    angular_resolution:
+        Angular resolution contribution.
+    wavelength_resolution:
+        Wavelength resolution contribution.
+    sample_size_resolution:
+        Sample size resolution contribution.
+
+    Returns
+    -------
+    :
+        Q resolution function.
+    """
+    return sc.sqrt(
+        (angular_resolution**2 + wavelength_resolution**2 + sample_size_resolution**2)
+        * Q**2
+    )
diff --git a/src/ess/estia/types.py b/src/ess/estia/types.py
new file mode 100644
index 00000000..024960f2
--- /dev/null
+++ b/src/ess/estia/types.py
@@ -0,0 +1,7 @@
+from typing import NewType
+
+import scipp as sc
+
+WavelengthResolution = NewType("WavelengthResolution", sc.Variable)
+AngularResolution = NewType("AngularResolution", sc.Variable)
+SampleSizeResolution = NewType("SampleSizeResolution", sc.Variable)
diff --git a/src/ess/estia/workflow.py b/src/ess/estia/workflow.py
new file mode 100644
index 00000000..8b137891
--- /dev/null
+++ b/src/ess/estia/workflow.py
@@ -0,0 +1 @@
+