Add involute surface to support HFIR simulation

[VHLM2001]

Adds involute surface to Orange, with additional solver and crests.

Involutes are curves are created by unwinding a chord from a parent shape.

The involute implemented here is constructed from a circle and can be made to be clockwise (negative) or anti-clockwise (positive).

This involute is the same type of that found in HFIR, and is necessary for building accurate models.

Involutes can be represented by the following equation:

$$\vec{r} = r_b\begin{bmatrix}\cos\left(t+a\right)+t\sin\left(t+a\right) \\ \sin\left(t+a\right)-t\cos\left(t+a\right) \\ 0 \end{bmatrix}$$

Where $\vec{r}$ is the parametric form of the involute, $r_b$ is the radius of the circle of involute, $a$ is the angular displacement of the involute from the x-axis, and $t$ is the angle from the displacement angle on describing the tangent that is used to construct the involute.

To determine the $a$ for this involute the tangent used for the construction of the involute for the point must be obtained. To obtain points instruction for the tangents that cross through that point, using the circle which the involute is constructed from and a second circle constructed with the radius being the distance from the origin to the point center point is needed, and shown bellow.

While the parameters define the curve itself have no restriction, except for the radius of involute being positive, the involute needs to be bounded for there to be finite solutions to the roots of the intersection points. Additionally this bound cannot exceed an interval of size of 2$\pi$ to be able to determine if whether a particle is in, out, or on the involute surface. Lastly the involute geometry is fixed to be z-aligned.

To determine the sense of a point, first it is determined if it is within the radial conversion of the bounds given by:

$$r=\sqrt{t^2+1}$$

Then it is checked if the particle is on the surface. After which it is possible to check if a particle is inside the region by finding if there exists an $a$ in the interval of $[a_0, a_0+t_{max}]$, where $a_0$ is the displacement angle for the reference involute and $t_{max}$ is the upper bound, for which an involute that passes through the point can be constructed.

This tangent can be obtained with circle of involute and constructing a second circle centered at the point with a radius equal to points distance from the origin shown bellow.

Then to find the $x\prime$ coordinate of the intersection of the intersection the following relation is used:

$$ x\prime = \frac{d^2-r^2+r_b^2}{2d} $$

Since $r=d$ the relation is simplified to:

$$ x\prime = \frac{r_b^2}{2r} $$

Then the $y\prime$ coordinate can be obtained from (positive for positive involutes and negative for negative involutes):

$$ y\prime = \pm \sqrt{r_b^2-{x\prime}^2} $$

Then the coordinates can be transformed to the $(x,y)$ by rotating the coordinates by the angle between the $x\prime$ -axis and the $x$-axis.

To calculate the intersections a root function of a line and involute is obtained bellow:

$$f(t) = r_b \left( v\cos\left(t+a\right) - u\sin\left(t+a\right) + t\left[\sin\left(t+a\right) + u\cos\left(t+a\right) \right] \right)-xu+yu$$

Where $u$ is the $x$ component of the direction of the particle, $v$ is the $y$ component of the direction of the particle, $x$ is the $x$ coordinate of the particle, and $y$ is the $y$ coordinate of the particle. Since this equation has no analytical solution an iteration scheme must be used. The method applied here is the Regular False Iteration due to its stability and described by:

$$ c_k= \frac{a_kf(b_k)-b_kf(a_k)}{f(b_k)-f(a_k) }$$

$c_k$ is used to update the bounds and find the root throughout the iteration, while initial values of $a_k$ and $b_k$ are the bounds the upper and lower bounds of the roots give by:

$$ \text{for } f(t)=0; t \in \{0, \beta - a, \beta - a - \pi, \beta - a + \pi, \beta - a -2\pi, \beta -a + 2\pi, ...\}$$

Where,

$$ \beta = \arctan\left(\frac{-v}{u}\right) $$

Lastly to calculate the outward normal of a point on the involute the following equation can be used:

$$ \hat{n} = s \begin{bmatrix} \sin\left(t+a\right) \\ -\cos\left(t+a\right) \\ 0 \end{bmatrix} $$

Where $s$ is the sign of the involute.

[sethrj]

Thanks for the contribution @VHLM2001 ! I'll be up at ORNL today and we can do a very rough first pass overview. In the meantime, can you please update the description to include:

• What is an involute and why does it matter than we can track on it?
• Are there any assumptions or restrictions about the surface?
• What methods are used to perform the intersection?

[sethrj]

@VHLM2001 can you please:

• Run git-clang-format $(git merge-base HEAD upstream/develop) and commit the result
• Run scripts/dev/install-commit-hooks.sh to install an auto-formatter for future commits
• Build with -Werror -Wall -Wextra -pedantic
• Merge upstream/develop

Thanks!

[tarapandya]

@VHLM2001 It would be good to include the equations for the involute in the description of the PR as well as a reference to the specific Regular Falsi iteration scheme you are using.

[VHLM2001]

Add involute images

[tarapandya]

@VHLM2001 Thanks for these updates. I looked through your tests. Be sure to reference the python files you added in the test files so someone knows they exists and what they are for. Also, looks like you have a lot of spelling typos in your test comments, please fix.

@sethrj Looks like this PR is ready for a more thorough review

[sethrj]

@VHLM2001 thanks for posting those plots! Would it be possible to regenerate them by calling ax.set_aspect("equal")? (so circles look like circles). I'll start taking another look at this today. Thanks!

[sethrj]

@VHLM2001 I've updated your branch through github by merging the latest upstream/develop into it. Please run git pull --rebase on your local machine. Thanks!

[sethrj]

@VHLM2001 And also take a look at and fix the CI failures: https://github.com/celeritas-project/celeritas/actions/runs/9810411562/job/27090497812?pr=1295 (or the "actions" tab at the top of this page). Whenever you push, those will re-run, and it may take a couple of iterations until all the tests pass, since there are several combinations of configurations.