Build a Custom Diffractometer Geometry#

This guide walks through defining a diffractometer geometry that is not one of the 10 pre-built presets. Each design decision is explained as it arises, with notes on what happens if you get it wrong.

The worked example builds a three-circle powder diffractometer with a theta-2theta arm and a sample spinner (phi). This is simple enough to follow but exercises every decision: basis, axis signs, parent chains, roles, and modes.

Tip

If your instrument matches one of the pre-built geometries (psic(), fourcv(), etc.), use the factory function instead. This guide is for instruments that differ from the presets.


The instrument#

Imagine a laboratory powder diffractometer with three rotary stages:

  1. theta — rotates the sample about the vertical axis (right-handed: positive angles rotate the sample face toward the beam, conventional laboratory sense).

  2. phi — a sample spinner mounted on top of theta, rotating about the axis perpendicular to the sample face (which, when theta = 0, points along the beam). Right-handed.

  3. ttheta — the detector arm, rotating about the same vertical axis as theta but mechanically independent. Right-handed.

The scattering plane is horizontal (both theta and ttheta rotate about the vertical axis, so the scattered beam sweeps through the horizontal plane).


Step 1 — Choose a basis#

The basis maps the three physical directions to Cartesian axes. Any right-handed orthogonal mapping works (see Choose and Understand Basis Vectors).

For a laboratory instrument the Busing & Levy convention is natural:

import numpy as np
import ad_hoc_diffractometer as ahd
from ad_hoc_diffractometer.factories import BASIS_BL
from ad_hoc_diffractometer.stage import Stage

basis = BASIS_BL
VERTICAL     = basis["vertical"]      # +z  [0, 0, 1]
LONGITUDINAL = basis["longitudinal"]  # +y  [0, 1, 0]
TRANSVERSE   = basis["transverse"]    # +x  [1, 0, 0]

Note

The basis choice does not affect the physics — only the numerical representation of axis vectors and matrices. Choosing a different basis produces the same motor angles from forward(). See Choose and Understand Basis Vectors for details.


Step 2 — Determine axis vectors from physical observation#

Stand at the instrument and observe each stage:

  1. What physical direction is the rotation axis? Look at the rotation axis of the stage. Theta and ttheta both rotate about the vertical axis. Phi rotates about the longitudinal axis (along the beam) when theta = 0.

  2. What is the handedness? Curl the fingers of your right hand around the axis so they point in the direction of positive rotation. If your thumb points along the positive physical direction, the rotation is right-handed (+). If it points opposite, the rotation is left-handed (-).

For this instrument:

Stage

Physical axis

Handedness

Axis vector

theta

vertical

right-handed

+VERTICAL

phi

longitudinal

right-handed

+LONGITUDINAL

ttheta

vertical

right-handed

+VERTICAL

Important

Getting the sign wrong inverts the rotation sense. The symptom is that forward() returns angles with the wrong sign — e.g. theta = -20 instead of +20 for a reflection you know requires positive rotation. If you see this, negate the axis vector for that stage.


Step 3 — Define parent relationships#

The parent argument on each Stage says which stage this one physically sits on:

  • parent=None — the stage is mounted on the fixed laboratory frame (the floor, the optical table, etc.).

  • parent="theta" — the stage is mounted on top of the theta stage; when theta rotates, this stage rotates with it.

For this instrument:

  • theta sits on the floor: parent=None

  • phi sits on top of theta: parent="theta"

  • ttheta is the detector arm, mechanically independent of the sample stages: parent=None

Important

Parent relationships encode the mechanical coupling between stages. If you set parent=None on a stage that is actually mounted on another stage, the forward solver will treat them as independent and compute incorrect angles.

Conversely, if two stacks share a common base stage (as in the zaxis() geometry, where the sample and detector both sit on the alpha stage), then both stacks must list that shared stage as their parent.


Step 4 — Assign roles#

Each stage has a role string:

  • "sample" — stages that orient the sample. Their rotation matrices combine to form the sample rotation matrix Z.

  • "detector" — stages that position the detector. Their rotation matrices combine to form the detector rotation matrix D.

Other roles (e.g. "analyzer", "polarizer") are accepted but are not used by the forward/inverse solvers.

For this instrument:

Stage

Role

theta

"sample"

phi

"sample"

ttheta

"detector"

Important

Assigning the wrong role causes the stage to be placed in the wrong rotation matrix (Z vs D). The symptom is that forward() produces solutions that do not satisfy the Bragg condition — inverse() of the returned angles does not give back the original (h, k, l).


Step 5 — Create the stages#

Putting the decisions together:

stages = [
    Stage("theta",  +VERTICAL,     parent=None,     role="sample"),
    Stage("phi",    +LONGITUDINAL, parent="theta",  role="sample"),
    Stage("ttheta", +VERTICAL,     parent=None,     role="detector"),
]

The order of stages in the list does not matter — the package builds the stacking order from the parent references. By convention, list sample stages first (floor-most first), then detector stages.


Step 6 — Define diffraction modes#

A ConstraintSet resolves the free degrees of freedom left after the Bragg condition is satisfied. An N-stage geometry has N - 3 free DOF:

  • 3-circle: 3 - 3 = 0 free DOF (no constraints needed)

  • 4-circle: 4 - 3 = 1 constraint

  • 6-circle: 6 - 3 = 3 constraints

This 3-circle geometry has zero free DOF, so each mode needs zero constraints. However, modes are still useful for documenting the intended operating condition and for the computed field (which lists the stages the solver writes to):

modes = {
    "theta_ttheta": ahd.ConstraintSet(
        [],
        computed=["theta", "phi", "ttheta"],
    ),
}

For geometries with free DOF, you must add constraints. See Work with Constraints and Diffraction Modes for the full constraint framework. Common patterns:

  • Bisecting: use BisectConstraint to tie a sample stage to half the detector angle (e.g. theta = ttheta / 2).

  • Fixed angle: use SampleConstraint to hold a stage at a declared value (e.g. phi = 0).

  • Frozen detector: use DetectorConstraint to hold a detector stage fixed (e.g. nu = 0 in psic).

Tip

The number of constraints must equal N - 3. Too few and the solver is under-determined; too many and it is over-determined. The package does not enforce this at construction time — it checks at solve time via is_fully_constrained().


Step 7 — Assemble the geometry#

Pass everything to AdHocDiffractometer:

g = ahd.AdHocDiffractometer(
    name="powder3c",
    stages=stages,
    basis=basis,
    description="Three-circle powder diffractometer, horizontal scattering plane",
    modes=modes,
    default_mode="theta_ttheta",
)

The constructor validates:

  • Stage names are unique

  • Parent references point to existing stages

  • No cycles in the parent graph

  • The basis contains exactly 3 mutually orthogonal, non-zero vectors

If any check fails, a ValueError is raised with a descriptive message.


Step 8 — Verify the geometry#

g.wavelength = 1.5406  # Cu K-alpha, Angstrom
g.sample.lattice = ahd.Lattice(a=5.431)  # cubic silicon
ahd.ub_identity(g.sample)

# Inspect
g.summary()

# Forward calculation: find motor angles for (1, 1, 1)
solutions = g.forward(1, 1, 1)
for sol in solutions:
    print(sol)

# Round-trip check: inverse should recover (1, 1, 1)
for sol in solutions:
    hkl = g.inverse(sol)
    print(f"  -> hkl = ({hkl[0]:.4f}, {hkl[1]:.4f}, {hkl[2]:.4f})")

If the round-trip does not recover the original (h, k, l) to within floating-point precision, check:

  1. Are the axis signs correct? (Step 2)

  2. Are the parent relationships correct? (Step 3)

  3. Are the roles correct? (Step 4)


Optional: register as a factory function#

To make the geometry available via list_geometries() and make_geometry(), wrap it in a function decorated with register_geometry():

from ad_hoc_diffractometer import register_geometry
from ad_hoc_diffractometer.factories import BASIS_BL

@register_geometry
def powder3c(basis=BASIS_BL):
    VERTICAL     = basis["vertical"]
    LONGITUDINAL = basis["longitudinal"]

    stages = [
        Stage("theta",  +VERTICAL,     parent=None,    role="sample"),
        Stage("phi",    +LONGITUDINAL, parent="theta", role="sample"),
        Stage("ttheta", +VERTICAL,     parent=None,    role="detector"),
    ]
    modes = {
        "theta_ttheta": ahd.ConstraintSet(
            [],
            computed=["theta", "phi", "ttheta"],
        ),
    }
    return ahd.AdHocDiffractometer(
        name="powder3c",
        stages=stages,
        basis=basis,
        description="Three-circle powder diffractometer",
        modes=modes,
        default_mode="theta_ttheta",
    )

The function name becomes the geometry name in the registry. Names must be unique — duplicates raise ValueError.


Advanced: shared base stages#

Some geometries have sample and detector stages that share a common base (e.g. the zaxis() geometry, where both the sample Z stage and the detector delta stage sit on the alpha stage). In this pattern:

stages = [
    Stage("alpha", +VERTICAL,     parent=None,    role="sample"),
    Stage("Z",     +LONGITUDINAL, parent="alpha", role="sample"),
    Stage("delta", -TRANSVERSE,   parent="alpha", role="detector"),
    Stage("gamma", +VERTICAL,     parent="delta", role="detector"),
]

The alpha stage is listed with role="sample" because it contributes to the sample rotation matrix Z. The detector stages (delta, gamma) reference parent="alpha" to express the mechanical coupling: when alpha rotates, the detector arm moves with it.

This pattern is used for surface diffraction geometries where the incidence angle is controlled by a shared tilt stage.


Decision checklist#

Use this checklist when defining a new geometry:

  • [ ] Basis — which physical direction maps to which Cartesian axis? (See Choose and Understand Basis Vectors)

  • [ ] Axis vector for each stage — what physical direction is the rotation axis, and is the rotation right-handed (+) or left-handed (-)?

  • [ ] Parent chain — which stage sits on which? Are sample and detector stacks independent or do they share a base stage?

  • [ ] Roles — is each stage "sample" or "detector"?

  • [ ] Modes — how many constraints are needed (N - 3)? Which constraint types resolve the free DOF?

  • [ ] Verification — does inverse(forward(h, k, l)) recover (h, k, l)?


See also#