SpECTRE
v2023.10.11

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 > ¢er_one, const std::array< double, 3 > ¢er_two, double radius_one, double radius_two, double z_plane_one, double z_plane_two)  
UniformCylindricalEndcap (UniformCylindricalEndcap &&)=default  
UniformCylindricalEndcap (const UniformCylindricalEndcap &)=default  
UniformCylindricalEndcap &  operator= (const UniformCylindricalEndcap &)=default 
UniformCylindricalEndcap &  operator= (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::NoFrame >  jacobian (const std::array< T, 3 > &source_coords) const 
template<typename T >  
tnsr::Ij< tt::remove_cvref_wrap_t< T >, 3, Frame::NoFrame >  inv_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) 
Map from 3D unit right cylinder to a volume that connects portions of two spherical surfaces.
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((zC_1^2)/R_1)\) is uniform in \(\bar{\rho} = \sqrt{\bar{x}^2+\bar{y}^2}\) and the angular coordinate \(\phi_1 = \atan((yC_1^1)/(xC_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((zC_2^2)/R_2)\) is uniform in \(\bar\rho\) and the angular coordinate \(\phi_2 = \atan((yC_2^1)/(xC_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()
maps \((\bar{x},\bar{y},\bar{z})\) to \((x,y,z)\) according to
\begin{align} x^0 &= C_1^0+\lambda(C_2^0C_1^0) + \cos\phi\left(R_1\sin\theta_1 + \lambda(R_2\sin\theta_2R_1\sin\theta_1)\right), \label{eq:x0} \\ x^1 &= C_1^1+\lambda(C_2^1C_1^1) + \sin\phi\left(R_1\sin\theta_1 + \lambda(R_2\sin\theta_2R_1\sin\theta_1)\right), \label{eq:x1} \\ x^2 &= C_1^2+\lambda(C_2^2C_1^2) + R_1\cos\theta_1 + \lambda(R_2\cos\theta_2R_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.
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^2C_1^2 + R_2\cos\theta_2R_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^0C_1^0\lambda(C_2^0C_1^0)\right)^2+ \left(x^1C_1^1\lambda(C_2^1C_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 onedimensional rootfinding 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^1C_1^1\lambda(C_2^1C_1^1)}{x^0C_1^0\lambda(C_2^0C_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}
We solve \(Q(\bar{\rho})=0\) numerically for \(\bar{\rho}\), where \(Q(\bar{\rho})\) is given by Eq. ( \(\ref{eq:defQ}\)).
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^2C_1^2 \geq R_1, \\ \displaystyle \frac{1}{\theta_{1 \mathrm{max}}} \cos^{1}\left(\frac{x^2C_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^2C_2^2 \leq R_2\cos\theta_{2 \mathrm{max}}, \\ \displaystyle \frac{1}{\theta_{2 \mathrm{max}}} \cos^{1}\left(\frac{x^2C_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}}\).
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 NewtonRaphson 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^0C_1^0\lambda(C_2^0C_1^0)\right)(C_2^0C_1^0)+ \left(x^1C_1^1\lambda(C_2^1C_1^1)\right)(C_2^1C_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_2R_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^2C_1^2 + R_2\cos\theta_2R_1\cos\theta_1} \label{eq:dlambda_drhobar}. \end{align}
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 rootfinders (e.g. NewtonRaphson) 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 firstorder 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\).
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 defaultconstructed 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 defaultconstructed 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^0C_1^0) + \cos\varphi (r_1 + \tilde{\lambda} (r_2 r_1)),\\ y_c &= C_1^1 + \tilde{\lambda}(C_2^1C_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^0C_1^0),\\ \tilde{y} &= y  C_1^1  \tilde{\lambda} (C_2^1C_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_2r_1)\tilde{\lambda}. \end{align}
The inverse map can therefore reject any points that do not satisfy this criterion.
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 wellbehaved 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^0C_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^1C_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^2C_1^2 + R_2\cos\theta_2  R_1\cos\theta_1 \right). \end{align}
The inverse Jacobian is computed by numerically inverting the Jacobian.
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
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
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
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 \((xC_1^0)/(yC_1^1) = (C_1^0C_2^0)/(C_1^1C_2^1)\), then the map is singlevalued everywhere. We cannot numerically check every point in this oneparameter 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.