SpECTRE  v2024.03.19
domain::CoordinateMaps::UniformCylindricalEndcap Class Reference

Map from 3D unit right cylinder to a volume that connects portions of two spherical surfaces. More...

#include <UniformCylindricalEndcap.hpp>

Public Member Functions

 UniformCylindricalEndcap (const std::array< double, 3 > &center_one, const std::array< double, 3 > &center_two, double radius_one, double radius_two, double z_plane_one, double z_plane_two)
 
 UniformCylindricalEndcap (UniformCylindricalEndcap &&)=default
 
 UniformCylindricalEndcap (const UniformCylindricalEndcap &)=default
 
UniformCylindricalEndcapoperator= (const UniformCylindricalEndcap &)=default
 
UniformCylindricalEndcapoperator= (UniformCylindricalEndcap &&)=default
 
template<typename T >
std::array< tt::remove_cvref_wrap_t< T >, 3 > operator() (const std::array< T, 3 > &source_coords) const
 
std::optional< std::array< double, 3 > > inverse (const std::array< double, 3 > &target_coords) const
 The inverse function is only callable with doubles because the inverse might fail if called for a point out of range, and it is unclear what should happen if the inverse were to succeed for some points in a DataVector but fail for other points.
 
template<typename T >
tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFramejacobian (const std::array< T, 3 > &source_coords) const
 
template<typename T >
tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrameinv_jacobian (const std::array< T, 3 > &source_coords) const
 
void pup (PUP::er &p)
 

Static Public Member Functions

static bool is_identity ()
 

Static Public Attributes

static constexpr size_t dim = 3
 

Friends

bool operator== (const UniformCylindricalEndcap &lhs, const UniformCylindricalEndcap &rhs)
 

Detailed Description

Map from 3D unit right cylinder to a volume that connects portions of two spherical surfaces.

A cylinder maps to the shaded region.

Details

Consider two spheres with centers \(C_1\) and \(C_2\), and radii \(R_1\) and \(R_2\). Let sphere 1 be intersected by a plane normal to the \(z\) axis and located at \(z = z_{\mathrm{P}1}\), and let sphere 2 be intersected by a plane normal to the \(z\) axis and located at \(z = z_{\mathrm{P}2}\).

UniformCylindricalEndcap maps a 3D unit right cylinder (with coordinates \((\bar{x},\bar{y},\bar{z})\) such that \(-1\leq\bar{z}\leq 1\) and \(\bar{x}^2+\bar{y}^2 \leq 1\)) to the shaded area in the figure above (with coordinates \((x,y,z)\)). The "bottom" of the cylinder \(\bar{z}=-1\) is mapped to the portion of sphere 1 that has \(z \geq z_{\mathrm{P}1}\), and on this portion of the sphere the angular coordinate \(\theta_1 = \acos((z-C_1^2)/R_1)\) is uniform in \(\bar{\rho} = \sqrt{\bar{x}^2+\bar{y}^2}\) and the angular coordinate \(\phi_1 = \atan((y-C_1^1)/(x-C_1^0))\) is the same as \(\phi = \atan(\bar{y}/\bar{x})\). Likewise, the "top" of the cylinder \(\bar{z}=+1\) is mapped to the portion of sphere 2 that has \(z \geq z_{\mathrm{P}2}\), and on this portion of the sphere the angular coordinate \(\theta_2 = \acos((z-C_2^2)/R_2)\) is uniform in \(\bar\rho\) and the angular coordinate \(\phi_2 = \atan((y-C_2^1)/(x-C_2^0))\) is the same as \(\phi\).

UniformCylindricalEndcap is different from CylindricalEndcap because of the distribution of points on the spheres, and because for UniformCylindricalEndcap the mapped portion of both Sphere 1 and Sphere 2 are bounded by planes of constant \(z\), whereas for CylindricalEndcap only one of the mapped portions is bounded by a plane (except for specially chosen map parameters). Note that UniformCylindricalEndcap can be used to construct maps that connect an arbitrary number of nested spheres; this is not possible for CylindricalEndcap for more than 3 nested spheres because of this asymmetry between CylindricalEndcap's two spherical surfaces.

UniformCylindricalEndcap is intended to be composed with Wedge<2> maps to construct a portion of a cylindrical domain for a binary system.

UniformCylindricalEndcap can be used to construct a domain that is similar to, but not identical to, the one described briefly in the Appendix of [32]. UniformCylindricalEndcap is used to construct the Blocks analogous to those labeled 'CA wedge', 'EA wedge', 'CB wedge', 'EE wedge', and 'EB wedge' in Figure 20 of that paper.

UniformCylindricalEndcap provides the following functions:

operator()

operator() maps \((\bar{x},\bar{y},\bar{z})\) to \((x,y,z)\) according to

\begin{align} x^0 &= C_1^0+\lambda(C_2^0-C_1^0) + \cos\phi\left(R_1\sin\theta_1 + \lambda(R_2\sin\theta_2-R_1\sin\theta_1)\right), \label{eq:x0} \\ x^1 &= C_1^1+\lambda(C_2^1-C_1^1) + \sin\phi\left(R_1\sin\theta_1 + \lambda(R_2\sin\theta_2-R_1\sin\theta_1)\right), \label{eq:x1} \\ x^2 &= C_1^2+\lambda(C_2^2-C_1^2) + R_1\cos\theta_1 + \lambda(R_2\cos\theta_2-R_1\cos\theta_1) \label{eq:x2}. \end{align}

Here

\begin{align} \lambda &= \frac{\bar{z}+1}{2},\label{eq:lambdafromzbar}\\ \theta_1 &= \bar{\rho} \theta_{1 \mathrm{max}},\label{eq:deftheta1}\\ \theta_2 &= \bar{\rho} \theta_{2 \mathrm{max}},\label{eq:deftheta2}\\ \phi &= \atan(\bar{y}/\bar{x})\label{eq:defphi}, \end{align}

where \(\theta_{1 \mathrm{max}}\) and \(\theta_{2 \mathrm{max}}\) are defined by

\begin{align} \cos(\theta_{1\mathrm{max}}) &= (z_{\mathrm{P}1}-C_1^2)/R_1,\\ \cos(\theta_{2\mathrm{max}}) &= (z_{\mathrm{P}2}-C_2^2)/R_2, \end{align}

and

\begin{align} \bar{\rho} &= \sqrt{\bar{x}^2+\bar{y}^2}/\bar{R} \label{eq:defrhobar}, \end{align}

where \(\bar{R}\) is the radius of the cylinder in barred coordinates, which is always unity.

inverse

Given \((x,y,z)\) we want to find \((\bar{x},\bar{y},\bar{z})\). From Eq. ( \(\ref{eq:x2}\)) we can write \(\lambda\) as a function of \(\bar\rho\):

\begin{align} \lambda &= \frac{x^2 - C_1^2 - R_1\cos\theta_1} {C_2^2-C_1^2 + R_2\cos\theta_2-R_1\cos\theta_1} \label{eq:lambda_from_rho}. \end{align}

Then by eliminating \(\phi\) from Eqs. ( \(\ref{eq:x0}\)) and ( \(\ref{eq:x1}\)) we find that \(\bar{\rho}\) is the solution of \(Q(\bar{\rho})=0\), where

\begin{align} Q(\bar{\rho}) &= \left(x^0-C_1^0-\lambda(C_2^0-C_1^0)\right)^2+ \left(x^1-C_1^1-\lambda(C_2^1-C_1^1)\right)^2- \left((1-\lambda)R_1\sin\theta_1 + \lambda R_2\sin\theta_2\right)^2.\label{eq:defQ} \end{align}

Here \(\lambda\), \(\theta_1\), and \(\theta_2\) are functions of \(\bar{\rho}\) through Eqs. ( \(\ref{eq:lambda_from_rho}\)), ( \(\ref{eq:deftheta1}\)), and ( \(\ref{eq:deftheta2}\)).

We solve \(Q(\bar{\rho})=0\) numerically; it is a one-dimensional root-finding problem.

Once we have determined \(\bar{\rho}\), we then obtain \(\lambda\) from Eq. ( \(\ref{eq:lambda_from_rho}\)), and we obtain \(\phi\) from

\begin{align} \tan\phi &= \frac{x^1-C_1^1-\lambda(C_2^1-C_1^1)}{x^0-C_1^0-\lambda(C_2^0-C_1^0)}. \end{align}

Then \(\bar{z}\) is obtained from Eq. ( \(\ref{eq:lambdafromzbar}\)) and \(\bar{x}\) and \(\bar{y}\) are obtained from

\begin{align} \bar{x} &= \bar{\rho}\bar{R}\cos\phi,\\ \bar{y} &= \bar{\rho}\bar{R}\sin\phi. \end{align}

Considerations when root-finding.

We solve \(Q(\bar{\rho})=0\) numerically for \(\bar{\rho}\), where \(Q(\bar{\rho})\) is given by Eq. ( \(\ref{eq:defQ}\)).

min/max values of \(\bar{\rho}\):

Note that the root we care about must have \(0\leq\lambda\leq 1\); therefore from Eq. ( \(\ref{eq:lambda_from_rho}\)) we have

\begin{align} \bar{\rho}_{\mathrm{min}} &= \left\{\begin{array}{ll} 0 & \text{for } x^2-C_1^2 \geq R_1, \\ \displaystyle \frac{1}{\theta_{1 \mathrm{max}}} \cos^{-1}\left(\frac{x^2-C_1^2}{R_1}\right) & \text{otherwise} \end{array}\right.\label{eq:rhobarmin}\\ \bar{\rho}_{\mathrm{max}} &= \left\{\begin{array}{ll} 1 & \text{for } x^2-C_2^2 \leq R_2\cos\theta_{2 \mathrm{max}}, \\ \displaystyle \frac{1}{\theta_{2 \mathrm{max}}} \cos^{-1}\left(\frac{x^2-C_2^2}{R_2}\right) & \text{otherwise} \end{array}\right.\label{eq:rhobarmax} \end{align}

so we look for a root only between \(\bar{\rho}_{\mathrm{min}}\) and \(\bar{\rho}_{\mathrm{max}}\).

Roots within roundoff of endpoints:

Sometimes a root is within roundoff of \(\bar{\rho}_{\mathrm{min}}\) or \(\bar{\rho}_{\mathrm{max}}\). This tends to happen at points on the boundary of the mapped region. In this case, the root might not be bracketed by \([\bar{\rho}_{\mathrm{min}},\bar{\rho}_{\mathrm{max}}]\) if the root is slightly outside that interval. If we find that \(Q(\bar{\rho}_{\mathrm{min}})\) is near zero but has the wrong sign, then we slightly expand the interval as follows:

\begin{align} \bar{\rho}_{\mathrm{min}} \to \bar{\rho}_{\mathrm{min}} - 2 \frac{Q(\bar{\rho}_{\mathrm{min}})}{Q'(\bar{\rho}_{\mathrm{min}})}, \end{align}

where \(Q'(\bar{\rho}_{\mathrm{min}})\) is the derivative of the function in Eq. ( \(\ref{eq:defQ}\)). Note that without the factor of 2, this is a Newton-Raphson step; the factor of 2 is there to overcompensate so that the new \(\bar{\rho}_{\mathrm{min}}\) brackets the root.

Similarly, if it is found that \(Q(\bar{\rho}_{\mathrm{max}})\) is near zero but has the wrong sign so that the root is not bracketed, then the same formula is used to expand the interval near \(\bar{\rho} = \bar{\rho}_{\mathrm{max}}\) to bracket the root.

Note that by differentiating Eqs. ( \(\ref{eq:defQ}\)) and ( \(\ref{eq:lambda_from_rho}\)), one obtains

\begin{align} Q'(\bar{\rho}) =& -2 \frac{d\lambda}{d\bar{\rho}}\left[ \left(x^0-C_1^0-\lambda(C_2^0-C_1^0)\right)(C_2^0-C_1^0)+ \left(x^1-C_1^1-\lambda(C_2^1-C_1^1)\right)(C_2^1-C_1^1) \right]\nonumber \\ & -2 \left((1-\lambda)R_1\sin\theta_1 + \lambda R_2\sin\theta_2\right) \left[ \frac{d\lambda}{d\bar{\rho}} (R_2\sin\theta_2-R_1\sin\theta_1) +(1-\lambda)R_1\theta_{1 \mathrm{max}}\cos\theta_1 +\lambda R_2\theta_{2 \mathrm{max}}\cos\theta_2 \right], \label{eq:defQderiv} \end{align}

where

\begin{align} \frac{d\lambda}{d\bar{\rho}} &= \frac{(1-\lambda)R_1\theta_{1 \mathrm{max}}\sin\theta_1 +\lambda R_2\theta_{2 \mathrm{max}}\sin\theta_2} {C_2^2-C_1^2 + R_2\cos\theta_2-R_1\cos\theta_1} \label{eq:dlambda_drhobar}. \end{align}

Roots within roundoff of \(\bar{\rho}=0\) or \(\bar{\rho}=1\):

For some points on the boundary of the mapped domain, the root will be within roundoff of \(\bar{\rho}=0\) or \(\bar{\rho}=1\). Here it does not always make sense to expand the range of the map if the root fails (by roundoff) to be bracketed, as is done above. Furthermore, when \(\bar{\rho}=0\) is a root it turns out that both \(Q(\bar{\rho})=0\) and \(Q'(\bar{\rho})=0\) for \(\bar{\rho}=0\), so some root-finders (e.g. Newton-Raphson) have difficulty converging. Therefore the cases where the root is within roundoff of \(\bar{\rho}=0\) or \(\bar{\rho}=1\) are treated separately.

These cases are detected by comparing terms in the first-order power series of \(Q(\bar{\rho})=0\) when expanded about \(\bar{\rho}=0\) or \(\bar{\rho}=1\). When one of these cases is recognized, the root is returned as either \(\bar{\rho}=0\) or \(\bar{\rho}=1\).

Quick rejection of points out of range of the map.

It is expected that inverse() will often be passed points \((x,y,z)\) that are out of the range of the map; in this case inverse() returns a default-constructed std::optional<std::array<double, 3>>. To avoid the difficulty and expense of attempting to solve \(Q(\bar{\rho})=0\) numerically for such points (and then having this solution fail), it is useful to quickly reject points \((x,y,z)\) that are outside the range of the map.

Any point in the range of the map must be inside or on sphere 2, and it must be outside or on sphere 1, so the inverse map can immediately return a default-constructed std::optional<std::array<double, 3>> for a point that does not satisfy these conditions.

Likewise, the inverse map can immediately reject any point with \(z < z_{\mathrm{P}1}\).

Finally, consider the circle \(S_1\) defining the intersection of sphere 1 and the plane \(z = z_{\mathrm{P}1}\); this circle has radius \(r_1 = R_1 \sin\theta_{1\mathrm{max}}\). Similarly, the circle \(S_2\) defining the intersection of sphere 2 and the plane \(z = z_{\mathrm{P}2}\) has radius \(r_2 = R_2 \sin\theta_{2\mathrm{max}}\). Now consider the cone that passes through these two circles. A point in the range of the map must be inside or on this cone. The cone can be defined parametrically as

\begin{align} x_c &= C_1^0 + \tilde{\lambda}(C_2^0-C_1^0) + \cos\varphi (r_1 + \tilde{\lambda} (r_2 -r_1)),\\ y_c &= C_1^1 + \tilde{\lambda}(C_2^1-C_1^1),+ \sin\varphi (r_1 + \tilde{\lambda} (r_2 -r_1)),\\ z_c &= C_1^2 + R_1 \cos\theta_{1\mathrm{max}} + \tilde{\lambda}(C_2^2 + R_2 \cos\theta_{2\mathrm{max}} - C_1^2 - R_1 \cos\theta_{1\mathrm{max}}), \end{align}

where \((x_c,y_c,z_c)\) is a point on the cone, and the two parameters defining a point on the cone are the angle \(\varphi\) around the cone and the parameter \(\tilde{\lambda}\), which is defined to be zero on \(S_1\) and unity on \(S_2\).

Given an arbitrary point \((x, y, z)\), we can determine whether or not that point is inside the cone as follows. First determine

\begin{align} \tilde{\lambda} &= \frac{z - C_1^2 - R_1 \cos\theta_{1\mathrm{max}}} {C_2^2+ R_2 \cos\theta_{2\mathrm{max}}- C_1^2- R_1 \cos\theta_{1\mathrm{max}}}, \\ \tilde{x} &= x - C_1^0 - \tilde{\lambda} (C_2^0-C_1^0),\\ \tilde{y} &= y - C_1^1 - \tilde{\lambda} (C_2^1-C_1^1).\\ \end{align}

Then the condition for the point to be inside or on the cone is

\begin{align} \sqrt{\tilde{x}^2+\tilde{y}^2} \le r_1 + (r_2-r_1)\tilde{\lambda}. \end{align}

The inverse map can therefore reject any points that do not satisfy this criterion.

jacobian

One can rewrite Eqs.( \(\ref{eq:x0}\)) through ( \(\ref{eq:x2}\)) as

\begin{align} x^0 &= \frac{1}{2}\left((1-\bar{z})C_1^0+ (1+\bar{z})C_2^0\right) + \frac{\bar{x}}{2}\left( (1-\bar{z}) R_1 S(\bar{\rho},\theta_{1 \mathrm{max}}) + (1+\bar{z}) R_2 S(\bar{\rho},\theta_{2 \mathrm{max}}) \right), \label{eq:x0alt} \\ x^1 &= \frac{1}{2}\left((1-\bar{z})C_1^1 + (1+\bar{z})C_2^1\right) + \frac{\bar{y}}{2}\left( (1-\bar{z})R_1 S(\bar{\rho},\theta_{1 \mathrm{max}}) + (1+\bar{z})R_2 S(\bar{\rho},\theta_{2 \mathrm{max}}) \right), \label{eq:x1alt} \\ x^2 &= \frac{1}{2}\left((1-\bar{z})C_1^2 + (1+\bar{z})C_2^2\right) + \left( (1-\bar{z})R_1 \cos\theta_1 + (1+\bar{z})R_2 \cos\theta_2 \right), \label{eq:x2alt} \\ \end{align}

where we have used Eq. ( \(\ref{eq:lambdafromzbar}\)) to eliminate \(\lambda\) in favor of \(\bar{z}\), and where we have defined the function

\begin{align} S(\bar{\rho},a) = \frac{\sin(\bar{\rho} a)}{\bar{\rho}}. \label{eq:Sdef} \end{align}

Note that \(S(\bar{\rho},a)\) is finite as \(\bar{\rho}\) approaches zero, and in the code we must take care that everything remains well-behaved in that limit.

Then differentiating Eqs. ( \(\ref{eq:x0alt}\)) and ( \(\ref{eq:x1alt}\)) with respect to \(\bar{x}\) and \(\bar{y}\), taking into account the dependence of \(\bar{\rho}\) on \(\bar{x}\) and \(\bar{y}\) from Eq. ( \(\ref{eq:defrhobar}\)), we find:

\begin{align} \frac{\partial x^0}{\partial \bar{x}} &= \frac{1}{2}\left( (1-\bar{z}) R_1 S(\bar{\rho},\theta_{1 \mathrm{max}}) + (1+\bar{z}) R_2 S(\bar{\rho},\theta_{2 \mathrm{max}}) \right) + \frac{\bar{x}^2}{2\bar{\rho}} \left[ (1-\bar{z}) R_1 S'(\bar{\rho},\theta_{1 \mathrm{max}}) + (1+\bar{z}) R_2 S'(\bar{\rho},\theta_{2 \mathrm{max}}) \right], \\ \frac{\partial x^1}{\partial \bar{y}} &= \frac{1}{2}\left( (1-\bar{z}) R_1 S(\bar{\rho},\theta_{1 \mathrm{max}}) + (1+\bar{z}) R_2 S(\bar{\rho},\theta_{2 \mathrm{max}}) \right) + \frac{\bar{y}^2}{2\bar{\rho}} \left[ (1-\bar{z}) R_1 S'(\bar{\rho},\theta_{1 \mathrm{max}}) + (1+\bar{z}) R_2 S'(\bar{\rho},\theta_{2 \mathrm{max}}) \right], \\ \frac{\partial x^0}{\partial \bar{y}} &= \frac{\bar{x}\bar{y}}{2\bar{\rho}}\left[ (1-\bar{z}) R_1 S'(\bar{\rho},\theta_{1 \mathrm{max}}) + (1+\bar{z}) R_2 S'(\bar{\rho},\theta_{2 \mathrm{max}}) \right], \\ \frac{\partial x^1}{\partial \bar{x}} &= \frac{\partial x^0}{\partial \bar{y}}, \end{align}

where \(S'(\bar{\rho},a)\) means the derivative of \(S(\bar{\rho},a)\) with respect to \(\bar\rho\). Note that \(S'(\bar{\rho},a)/\bar{\rho}\) approaches a constant value as \(\bar{\rho}\) approaches zero.

Differentiating Eq. ( \(\ref{eq:x2alt}\)) with respect to \(\bar{x}\) and \(\bar{y}\) we find

\begin{align} \frac{\partial x^2}{\partial \bar{x}} &= - \frac{\bar{x}}{2}\left[ (1-\bar{z}) R_1 \theta_{1 \mathrm{max}} S(\bar{\rho},\theta_{1 \mathrm{max}}) + (1+\bar{z}) R_2 \theta_{2 \mathrm{max}} S(\bar{\rho},\theta_{2 \mathrm{max}})- \right],\\ \frac{\partial x^2}{\partial \bar{y}} &= - \frac{\bar{y}}{2}\left[ (1-\bar{z}) R_1 \theta_{1 \mathrm{max}} S(\bar{\rho},\theta_{1 \mathrm{max}}) + (1+\bar{z}) R_2 \theta_{2 \mathrm{max}} S(\bar{\rho},\theta_{2 \mathrm{max}})- \right]. \end{align}

Differentiating Eqs. ( \(\ref{eq:x0alt}\)) through ( \(\ref{eq:x2alt}\)) with respect to \(\bar{z}\) yields

\begin{align} \frac{\partial x^0}{\partial \bar{z}} &= \frac{1}{2}\left[ C_2^0-C_1^0 + \bar{x}\left(R_2 S(\bar{\rho},\theta_{2 \mathrm{max}}) - R_1 S(\bar{\rho},\theta_{1 \mathrm{max}})\right) \right],\\ \frac{\partial x^1}{\partial \bar{z}} &= \frac{1}{2}\left[ C_2^1-C_1^1 + \bar{y}\left(R_2 S(\bar{\rho},\theta_{2 \mathrm{max}}) - R_1 S(\bar{\rho},\theta_{1 \mathrm{max}})\right) \right],\\ \frac{\partial x^2}{\partial \bar{z}} &= \frac{1}{2}\left( C_2^2-C_1^2 + R_2\cos\theta_2 - R_1\cos\theta_1 \right). \end{align}

inv_jacobian

The inverse Jacobian is computed by numerically inverting the Jacobian.

Restrictions on map parameters

We demand that Sphere 1 is fully contained inside Sphere 2, and that the two spheres have at least some small separation between them. In particular, we demand that

  • \(0.98 R_2 >= R_1 + |C_1-C_2|\)

where 0.98 is a safety factor. It is possible to construct a valid map without this assumption, but the assumption simplifies the code, and the expected use cases obey this restriction.

We also demand that \(z_{\mathrm{P}1} <= z_{\mathrm{P}2} -0.04 R_2\), and that the z planes in the above figures lie above the centers of the corresponding spheres and are not too close to the centers or edges of those spheres; specificially, we demand that

  • \( 0.075\pi < \theta_{1 \mathrm{max}} < 0.45\pi\)
  • \( 0.075\pi < \theta_{2 \mathrm{max}} < 0.45\pi\)

Here 0.075 and 0.45 are safety factors. These restrictions are not strictly necessary but are made for simplicity and to ensure the accuracy of the inverse map (the inverse map becomes less accurate if the map parameters are extreme).

Consider the line segment \(L\) that connects a point on the circle \(S_1\) (the circle formed by the intersection of sphere 1 and the plane \(z=z_{\mathrm{P}1}\)) with the center of the circle \(S_1\). Consider another line segment \(L'\) that connects the same point on the circle \(S_1\) with the corresponding point on the circle \(S_2\) (the circle formed by the intersection of sphere 2 and the plane \(z=z_{\mathrm{P}2}\)). Now consider the angle between \(L\) and \(L'\), as measured from the interior of sphere 1, and Let \(\alpha\) be the minimum value of this angle over the circle. \(\alpha\) is shown in the figure above. If \(\alpha < \theta_{1 \mathrm{max}}\), then the line segment \(L'\) intersects the mapped portion of sphere 1 twice, so the map is multiply valued. Similarly, if \(\alpha < \theta_{2 \mathrm{max}}\), then the line segment \(L'\) intersects the mapped portion of sphere 2 twice, and again the map is multiply valued. Therefore we demand that the map parameters are such that

  • \(\alpha > 1.1 \theta_{1 \mathrm{max}}\)
  • \(\alpha > 1.1 \theta_{2 \mathrm{max}}\)

where 1.1 is a safety factor.

The condition on \(\alpha\) is guaranteed to provide an invertible map if \(C_1^0=C_2^0\) and \(C_1^1=C_2^1\). However, for \(C_1^0 \neq C_2^0\) or \(C_1^1\neq C_2^1\), even if the \(\alpha\) condition is satisfied, it is possible for two lines of constant \((\bar{x},\bar{y})\) (each line has different values of \((\bar{x},\bar{y})\)) to pass through the same point \((x,y,z)\) if those lines are not coplanar. This condition is difficult to check analytically, so we check it numerically. We have found empirically that if \(Q(\bar{\rho})\) from Eq. ( \(\ref{eq:defQ}\)) has only a single root between \(\bar{\rho}_{\mathrm{min}}\) and \(\bar{\rho}_{\mathrm{max}}\) for all points \((x,y,z)\) on the surface of sphere 1 with \(z\geq z_{\mathrm{P}1}\) and with \((x-C_1^0)/(y-C_1^1) = (C_1^0-C_2^0)/(C_1^1-C_2^1)\), then the map is single-valued everywhere. We cannot numerically check every point in this one-parameter family of points, but we demand that this condition is satisfied for a reasonably large number of points in this family. This check is not very expensive since it is done only once, in the constructor.


The documentation for this class was generated from the following file: