Implement torus

[sethrj]

Existing implementations

OpenMC

https://github.com/openmc-dev/openmc/pull/1933, https://github.com/openmc-dev/openmc/pull/2053

MCNP/Lava

Elliptical or circular torus with axis alignment:

$$(z - \bar z)^2 / b^2 + \left( \sqrt{(x - \bar x)^2 + (y - \bar y)^2} - a\right)^2 /c^2 - 1 = 0$$

where $\bar x, \bar y, \bar z, a, b, c$ specify the ellipse

$$ \frac{s^2}{b^2} + \frac{(r-a)^2}{c^2} = 1 $$

rotated about the $s$ axis in the $(r, s)$ cylindrical coordinate system whose center is $(\bar x, \bar y, \bar z)$. For a Z torus, $s = z - \bar z$ and $r = \sqrt{(x - \bar x)^2 + (y - \bar y)^2}$

• $a$ is the distance from the centerline of the torus to the center of its ellipse
• $b$ is the ellipse axis along the main axis, $c$ along the minor
• MCNP supports "degenerate" torus: like an ellipsoid with a bump taken out of it ($0 < a < c$)
• MCNP deletes the torus if $|a| < c$
• Lava multiplies through by $c^2$
• UVW transformation: {{x,y,z}, {y, z, x}, {z, x, y}}

Normal:

    double dx = pos[0] - sp->xc;
    double dy = pos[1] - sp->yc;
    double dz = pos[2] - sp->zc;
    double tyz = 1.0 - (sp->a * 1.0 / sqrt(dy*dy + dz*dz));

    nrm[0] = sp->c2_b2 * dx;
    nrm[1] = dy * tyz;
    nrm[2] = dz * tyz;
    normalize_vector(nrm);

Geant4/VecGeom

Implementations:

• Geant4: G4Torus
• Vecgeom surface (new): TorusImpl
• Vecgeom volume (old): TorusImplementation2

Notes:

• Only Z axis
• Only circular
• Only centered
• Probably not degenerate??
• Intersection of +torus, -torus, and wedge
• VecGeom accelerates using bounding cylindrical segment
• VecGeom 2 moves first to bounding cylinder, then to torus
• Geant4 prohibits degenerate (and nearly degenerate!) torus: fails if $a < c + \epsilon$
• Geant4 uses Jenkins-Traub polynomial root solver

Variable	MCNP	Geant4
Swept center radius	a	fRtor
Radius along axis	b	fRmax
Radius perp to axis	c	fRmax
Inner radius	—	fRmin

Possible simplifications

• Centered on origin: x, y, z are zero
• Circular cross section: b = c
• Degenerate $a=0$: convert to sphere or ellipsoid
• $c < a$: delete shape

Steps

• [ ] Implement Ferrari-Cardano quartic polynomial solver in orange/surf/detail, including comprehensive unit tests: #2010
• [ ] Implement Algorithm 1010 using the same unit test cases
• [ ] Add an Toroid surface to orange/surf, following the example of #1295
• [ ] Hookup the Toroid surface to rest of the code, following the example of #1342
• [ ] Add a Geant4 support for torii within orange/g4org/SolidConverter.