|
SpECTRE
v2024.05.11
|
Map from 3D unit right cylindrical shell to a volume that connects portions of two spherical surfaces. More...
#include <UniformCylindricalSide.hpp>
Public Member Functions | |
| UniformCylindricalSide (const std::array< double, 3 > ¢er_one, const std::array< double, 3 > ¢er_two, double radius_one, double radius_two, double z_plane_plus_one, double z_plane_minus_one, double z_plane_plus_two, double z_plane_minus_two) | |
| UniformCylindricalSide (UniformCylindricalSide &&)=default | |
| UniformCylindricalSide (const UniformCylindricalSide &)=default | |
| UniformCylindricalSide & | operator= (const UniformCylindricalSide &)=default |
| UniformCylindricalSide & | operator= (UniformCylindricalSide &&)=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 UniformCylindricalSide &lhs, const UniformCylindricalSide &rhs) |
Map from 3D unit right cylindrical shell 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\). Sphere 1 is assumed to be contained inside Sphere 2. Let sphere 1 be intersected by two planes normal to the \(z\) axis and located at \(z = z^{\pm}_{\mathrm{P}1}\), and let sphere 2 be intersected by two planes normal to the \(z\) axis and located at \(z = z^{\pm}_{\mathrm{P}2}\). Here we assume that \(z^{-}_{\mathrm{P}2} \leq z^{-}_{\mathrm{P}1}< z^{+}_{\mathrm{P}1} \leq z^{+}_{\mathrm{P}2}\).
UniformCylindricalSide maps a 3D unit right cylindrical shell (with coordinates \((\bar{x},\bar{y},\bar{z})\) such that \(-1\leq\bar{z}\leq 1\) and \(1 \leq \bar{x}^2+\bar{y}^2 \leq 4\), where the values of 1 and 2 for the inner and outer cylindrical radii are arbitrary choices but are required by UniformCylindricalSide) to the shaded area in the figure above (with coordinates \((x,y,z)\)). The "inner surface" of the cylindrical shell \(\bar{x}^2+\bar{y}^2=1\) is mapped to the portion of sphere 1 that has \(z^{-}_{\mathrm{P}1} \leq z \leq z^{+}_{\mathrm{P}1} \), and on this portion of the sphere the cosine of the polar angular coordinate \(\cos\theta_1 =(z-C_1^z)/R_1\) is uniform in \(\bar{z}\), and the angular coordinate \(\phi_1 = \atan((y-C_1^y)/(x-C_1^x))\) is the same as \(\phi = \atan(\bar{y}/\bar{x})\). Likewise, the "outer surface" of the cylindrical shell \(\bar{x}^2+\bar{y}^2=4\) is mapped to the portion of sphere 2 that has \(z^{-}_{\mathrm{P}2} \leq z \leq z^{+}_{\mathrm{P}2} \), and on this portion of the sphere the cosine of the azimuthal angular coordinate \(\cos\theta_2 = (z-C_2^z)/R_2\) is uniform in \(\bar{z}\), and the angular coordinate \(\phi_2 = \atan((y-C_2^y)/(x-C_2^x))\) is the same as \(\phi\).
UniformCylindricalSide is different from CylindricalSide because of the distribution of points on the spheres, and because for UniformCylindricalSide the mapped portion of both Sphere 1 and Sphere 2 are bounded by planes of constant \(z\), whereas for CylindricalSide only one of the mapped portions is bounded by a plane (except for specially chosen map parameters). Note that UniformCylindricalSide can be used to construct maps that connect an arbitrary number of nested spheres; this is not possible for CylindricalSide for more than 3 nested spheres because of this asymmetry between CylindricalSide's two spherical surfaces.
Note that the entire region between Sphere 1 and Sphere 2 can be covered by a single cylindrical shell (mapped using UniformCylindricalSide) and two cylinders (each mapped by UniformCylindricalEndcap).
UniformCylindricalSide is intended to be composed with Wedge<2> maps to construct a portion of a cylindrical domain for a binary system.
UniformCylindricalSide can be used to construct a domain that is similar to, but not identical to, the one described briefly in the Appendix of [32]. UniformCylindricalSide is used to construct the Blocks analogous to those labeled 'CA cylinder', 'EA cylinder', 'CB cylinder', 'EE cylinder', and 'EB cylinder' in Figure 20 of that paper.
UniformCylindricalSide provides the following functions:
operator() maps \((\bar{x},\bar{y},\bar{z})\) to \((x,y,z)\) according to
\begin{align} x &= C_1^x+\lambda(C_2^x-C_1^x) + \cos\phi\left(R_1\sin\theta_1 + \lambda(R_2\sin\theta_2-R_1\sin\theta_1)\right), \label{eq:x0} \\ y &= C_1^y+\lambda(C_2^y-C_1^y) + \sin\phi\left(R_1\sin\theta_1 + \lambda(R_2\sin\theta_2-R_1\sin\theta_1)\right), \label{eq:x1} \\ z &= C_1^z+\lambda(C_2^z-C_1^z) + 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 &= \bar{\rho}-1,\label{eq:lambdafromrhobar}\\ \cos\theta_1 &= \cos\theta_{1 \mathrm{max}} + \left(\cos\theta_{1 \mathrm{min}}-\cos\theta_{1 \mathrm{max}}\right) \frac{\bar{z}+1}{2}\label{eq:deftheta1}\\ \cos\theta_2 &= \cos\theta_{2 \mathrm{max}} + \left(\cos\theta_{2 \mathrm{min}}- \cos\theta_{2 \mathrm{max}}\right) \frac{\bar{z}+1}{2}\label{eq:deftheta2}\\ \phi &= \atan(\bar{y}/\bar{x})\label{eq:defphi}, \end{align}
where \(\theta_{1 \mathrm{min}}\), \(\theta_{2 \mathrm{min}}\), \(\theta_{1 \mathrm{max}}\), and \(\theta_{2 \mathrm{max}}\) are defined by
\begin{align} \label{eq:deftheta1min} \cos(\theta_{1\mathrm{min}}) &= (z^{+}_{\mathrm{P}1}-C_1^z)/R_1,\\ \cos(\theta_{1\mathrm{max}}) &= (z^{-}_{\mathrm{P}1}-C_1^z)/R_1,\\ \cos(\theta_{2\mathrm{min}}) &= (z^{+}_{\mathrm{P}2}-C_2^z)/R_2,\\ \label{eq:deftheta2max} \cos(\theta_{2\mathrm{max}}) &= (z^{-}_{\mathrm{P}2}-C_2^z)/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 inner radius of the cylindrical shell in barred coordinates, which is always unity.
Note that \(\theta_{1\mathrm{min}}<\theta_{1\mathrm{max}}\) but \(\cos\theta_{1\mathrm{min}}>\cos\theta_{1\mathrm{max}}\) (and same for sphere 2).
Also note that Eqs. ( \(\ref{eq:deftheta1}\)) and ( \(\ref{eq:deftheta2}\)) can be simplified using Eqs. ( \(\ref{eq:deftheta1min}\)- \(\ref{eq:deftheta2max}\)):
\begin{align} R_1\cos\theta_1 &= z^{-}_{\mathrm{P}1}-C_1^z +(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1}) \frac{\bar{z}+1}{2}\label{eq:deftheta1alt}\\ R_2\cos\theta_2 &= z^{-}_{\mathrm{P}2}-C_2^z +(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2}) \frac{\bar{z}+1}{2}\label{eq:deftheta2alt}\\ \end{align}
Given \((x,y,z)\) we want to find \((\bar{x},\bar{y},\bar{z})\). From Eqs. ( \(\ref{eq:x2}\)), ( \(\ref{eq:deftheta1alt}\)), and ( \(\ref{eq:deftheta2alt}\)) we can write \(\bar{z}\) as a function of \(\lambda\):
\begin{align} \frac{1+\bar{z}}{2} &= \frac{z + \lambda (z^{-}_{\mathrm{P}1}-z^{-}_{\mathrm{P}2}) - z^{-}_{\mathrm{P}1}} {(1-\lambda)(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1}) + \lambda(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2})} \label{eq:zbar_from_lambda}, \end{align}
Note that the denominator of Eq. ( \(\ref{eq:zbar_from_lambda}\)) is always positive because \(0\leq\lambda\leq 1\), \(z^{+}_{\mathrm{P}1}>z^{-}_{\mathrm{P}1}\), and \(z^{+}_{\mathrm{P}2}>z^{-}_{\mathrm{P}2}\).
By eliminating \(\phi\) from Eqs. ( \(\ref{eq:x0}\)) and ( \(\ref{eq:x1}\)) we find that \(\lambda\) is the solution of \(Q(\lambda)=0\), where
\begin{align} Q(\lambda) &= \left(x-C_1^x-\lambda(C_2^x-C_1^x)\right)^2+ \left(y-C_1^y-\lambda(C_2^y-C_1^y)\right)^2- \left((1-\lambda)R_1\sin\theta_1 + \lambda R_2\sin\theta_2\right)^2.\label{eq:defQ} \end{align}
Here \(\theta_1\) and \(\theta_2\) are functions of \(\bar{z}\) through Eqs. ( \(\ref{eq:deftheta1alt}\)) and ( \(\ref{eq:deftheta2alt}\)), and \(\bar{z}\) is a function of \(\lambda\) through Eq. ( \(\ref{eq:zbar_from_lambda}\)).
We solve \(Q(\lambda)=0\) numerically; it is a one-dimensional root-finding problem.
Once we have determined \(\lambda\), we then obtain \(\bar{z}\) from Eq. ( \(\ref{eq:zbar_from_lambda}\)), and we obtain \(\phi\) from
\begin{align} \tan\phi &= \frac{y-C_1^y-\lambda(C_2^y-C_1^y)}{x-C_1^x-\lambda(C_2^x-C_1^x)}. \end{align}
Then \(\bar{\rho}\) is obtained from Eq. ( \(\ref{eq:lambdafromrhobar}\)) 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(\lambda)=0\) numerically for \(\lambda\), where \(Q(\lambda)\) is given by Eq. ( \(\ref{eq:defQ}\)).
Note that the root we care about must have \(-1\leq\bar{z}\leq 1\); therefore from Eq. ( \(\ref{eq:zbar_from_lambda}\)) we have
\begin{align} \lambda_{\mathrm{min}} &= \mathrm{max}\left\{0, \frac{z-z^{+}_{\mathrm{P}1}} {(z^{+}_{\mathrm{P}2}-z^{+}_{\mathrm{P}1})}, \frac{z^{-}_{\mathrm{P}1}-z} {(z^{-}_{\mathrm{P}1}-z^{-}_{\mathrm{P}2})} \right\}\label{eq:lambdamin} \end{align}
In the case where \(z^{+}_{\mathrm{P}2}=z^{+}_{\mathrm{P}1}\) we treat the middle term in Eq.( \(\ref{eq:lambdamin}\)) as zero since in that case \(z-z^{+}_{\mathrm{P}1}\) can never be positive for \(x^2\) in the range of the map, and for \(z=z^{+}_{\mathrm{P}2}=z^{+}_{\mathrm{P}1}\) it turns out that ( \(\ref{eq:zbar_from_lambda}\)) places no restriction on \(\lambda_{\mathrm{min}}\). For the same reason, if \(z^{-}_{\mathrm{P}2}=z^{-}_{\mathrm{P}1}\) we treat the last term in Eq.( \(\ref{eq:lambdamin}\)) as zero.
We look for a root only between \(\lambda_{\mathrm{min}}\) and \(\lambda_{\mathrm{max}}=1\).
Sometimes a root is within roundoff of \(\lambda_{\mathrm{min}}\). In this case, the root might not be bracketed by \([\lambda_{\mathrm{min}},\lambda_{\mathrm{max}}]\) if the root is slightly outside that interval by roundoff error. If we find that \(Q(\lambda_{\mathrm{min}})\) is near zero but has the wrong sign, then we slightly expand the interval as follows:
\begin{align} \lambda_{\mathrm{min}} \to \lambda_{\mathrm{min}} - 2 \frac{Q(\lambda_{\mathrm{min}})}{Q'(\lambda_{\mathrm{min}})}, \end{align}
where \(Q'(\lambda_{\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 \(\lambda_{\mathrm{min}}\) brackets the root.
Note that by differentiating Eqs. ( \(\ref{eq:defQ}\)) and ( \(\ref{eq:zbar_from_lambda}\)), one obtains
\begin{align} Q'(\lambda) =& -2 \left[ \left(x-C_1^x-\lambda(C_2^x-C_1^x)\right)(C_2^x-C_1^x)+ \left(y-C_1^y-\lambda(C_2^y-C_1^y)\right)(C_2^y-C_1^y) \right]\nonumber \\ & -\left[ 2(R_2\sin\theta_2-R_1\sin\theta_1) -(1-\lambda)\cot\theta_1 (z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1}) \frac{d\bar{z}}{d\lambda} \right. \nonumber \\ & \left.\qquad -\lambda \cot\theta_2 (z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2}) \frac{d\bar{z}}{d\lambda} \right] \left((1-\lambda)R_1\sin\theta_1 + \lambda R_2\sin\theta_2\right), \label{eq:defQderiv} \end{align}
where
\begin{align} \frac{d\bar{z}}{d\lambda} &= \frac{(1-\bar{z})(z^{-}_{\mathrm{P}1}-z^{-}_{\mathrm{P}2}) -(1+\bar{z})(z^{+}_{\mathrm{P}2}-z^{+}_{\mathrm{P}1})} {(1-\lambda)(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1}) + \lambda(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2})} \label{eq:dzbar_dlambda}. \end{align}
A root within roundoff of \(\lambda_{\mathrm{max}}\) is treated similarly.
For some points on the boundary of the mapped domain, \(\lambda_{\mathrm{min}}\) will be within roundoff of \(\lambda=1\). We check explicitly for this case, and we compute the root as exactly \(\lambda=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 std::nullopt. To avoid the difficulty and expense of attempting to solve \(Q(\lambda)=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}2}\) or \(z > z^{+}_{\mathrm{P}2}\).
Finally, for \(z^{+}_{\mathrm{P}2}\neq z^{+}_{\mathrm{P}1}\), 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{min}}\). 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{min}}\). Now consider the cone that passes through these two circles. A point in the range of the map must be outside (where "outside" means farther from the \(z\) axis) or on this cone. The cone can be defined parametrically as
\begin{align} x_c &= C_1^x + \tilde{\lambda}(C_2^x-C_1^x) + \cos\varphi (r_1 + \tilde{\lambda} (r_2 -r_1)),\\ y_c &= C_1^y + \tilde{\lambda}(C_2^y-C_1^y),+ \sin\varphi (r_1 + \tilde{\lambda} (r_2 -r_1)),\\ z_c &= z^{+}_{\mathrm{P}1} + \tilde{\lambda}(z^{+}_{\mathrm{P}2}-z^{+}_{\mathrm{P}1}), \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 - z^{+}_{\mathrm{P}1}} {z^{+}_{\mathrm{P}2}-z^{+}_{\mathrm{P}1}}, \\ \tilde{x} &= x - C_1^x - \tilde{\lambda} (C_2^x-C_1^x),\\ \tilde{y} &= y - C_1^y - \tilde{\lambda} (C_2^y-C_1^y).\\ \end{align}
Then the condition for the point to be outside or on the cone is
\begin{align} \sqrt{\tilde{x}^2+\tilde{y}^2} \ge r_1 + (r_2-r_1)\tilde{\lambda}. \end{align}
The inverse map therefore rejects any points that do not satisfy this criterion. The cone criterion makes sense only for points with \(z\geq z^{+}_{\mathrm{P}1}\).
For \(z^{-}_{\mathrm{P}2} \neq z^{-}_{\mathrm{P}1}\), a similar cone can be constructed for the southern hemisphere. That cone passes through the circle \(S^{-}_1\) defining the intersection of sphere 1 and the plane \(z = z^{-}_{\mathrm{P}1}\) and the circle \(S^{-}_2\) defining the intersection of sphere 2 and the plane \(z = z^{-}_{\mathrm{P}2}\). The inverse map rejects any point that is inside that cone as well, provided that the point has \(z\leq z^{-}_{\mathrm{P}1}\). For points with \(z > z^{-}_{\mathrm{P}1}\) checking the cone criterion does not make sense.
From Eqs. ( \(\ref{eq:deftheta1alt}\)) and ( \(\ref{eq:deftheta2alt}\)) we see that \(\theta_1\) and \(\theta_2\) depend on \(\bar{z}\) and are independent of \(\bar{x}\) and \(\bar{y}\), and that
\begin{align} \frac{\partial (R_1\cos\theta_1)}{\partial\bar{z}} &= \frac{1}{2}(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1}), \label{eq:dcostheta1} \\ \frac{\partial (R_2\cos\theta_2)}{\partial\bar{z}} &= \frac{1}{2}(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2}), \label{eq:dcostheta2} \\ \frac{\partial (R_1\sin\theta_1)}{\partial\bar{z}} &= -\frac{1}{2}\cot\theta_1 (z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1}), \label{eq:dsintheta1} \\ \frac{\partial (R_2\sin\theta_2)}{\partial\bar{z}} &= -\frac{1}{2}\cot\theta_2(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2}), \label{eq:dsintheta2} \end{align}
Also, from Eqs. ( \(\ref{eq:defphi}\)) and ( \(\ref{eq:defrhobar}\)) we have
\begin{align} \frac{\partial\cos\phi}{\partial\bar{x}} &= \frac{\bar{y}^2}{\bar{R}^3\bar{\rho}^3}, \label{eq:dcosphidxbar} \\ \frac{\partial\cos\phi}{\partial\bar{y}} &= -\frac{\bar{x}\bar{y}}{\bar{R}^3\bar{\rho}^3}, \label{eq:dcosphidybar} \\ \frac{\partial\sin\phi}{\partial\bar{x}} &= -\frac{\bar{x}\bar{y}}{\bar{R}^3\bar{\rho}^3}, \label{eq:dsinphidxbar} \\ \frac{\partial\sin\phi}{\partial\bar{y}} &= \frac{\bar{x}^2}{\bar{R}^3\bar{\rho}^3}, \label{eq:dsinphidybar} \end{align}
and we know that \(\phi\) is independent of \(\bar{z}\).
Finally, from Eqs. ( \(\ref{eq:defrhobar}\)) and ( \(\ref{eq:lambdafromrhobar}\)) we have
\begin{align} \frac{\partial\lambda}{\partial\bar{x}} &= \frac{\bar{x}}{\bar{R}^2\bar{\rho}}, \label{eq:dlambdadxbar} \\ \frac{\partial\lambda}{\partial\bar{y}} &= \frac{\bar{y}}{\bar{R}^2\bar{\rho}}, \label{eq:dlambdadybar} \end{align}
with no dependence on \(\bar{z}\).
Putting these results together yields
\begin{align} \frac{\partial x^0}{\partial \bar{x}} &= \frac{\bar{y}^2}{\bar{\rho}^3\bar{R}^3}R_1\sin\theta_1 + (R_2\sin\theta_2-R_1\sin\theta_1) \frac{\lambda \bar{R}^2\bar\rho^2+\bar{x}^2}{\bar\rho^3\bar{R}^3} + \frac{\bar{x}}{\bar\rho\bar{R}^2}(C_2^x-C_1^x),\\ \frac{\partial x^0}{\partial \bar{y}} &= \frac{\bar{x}\bar{y}}{\bar{\rho}^3\bar{R}^3} (R_2\sin\theta_2-2 R_1\sin\theta_1) + \frac{\bar{y}}{\bar\rho\bar{R}^2}(C_2^x-C_1^x),\\ \frac{\partial x^0}{\partial \bar{z}} &= -\frac{1}{2}\frac{\bar{x}}{\bar\rho\bar{R}}\left[ \cot\theta_1(1-\lambda)(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1})+ \cot\theta_2\lambda(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2})\right],\\ \frac{\partial x^1}{\partial \bar{x}} &= \frac{\bar{x}\bar{y}}{\bar{\rho}^3\bar{R}^3} (R_2\sin\theta_2-2 R_1\sin\theta_1) + \frac{\bar{x}}{\bar\rho\bar{R}^2}(C_2^y-C_1^y),\\ \frac{\partial x^1}{\partial \bar{y}} &= \frac{\bar{x}^2}{\bar{\rho}^3\bar{R}^3}R_1\sin\theta_1 + (R_2\sin\theta_2-R_1\sin\theta_1) \frac{\lambda \bar{R}^2\bar\rho^2+\bar{y}^2}{\bar\rho^3\bar{R}^3} + \frac{\bar{y}}{\bar\rho\bar{R}^2}(C_2^y-C_1^y),\\ \frac{\partial x^1}{\partial \bar{z}} &= -\frac{1}{2}\frac{\bar{y}}{\bar\rho\bar{R}}\left[ \cot\theta_1(1-\lambda)(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1})+ \cot\theta_2\lambda(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2})\right],\\ \frac{\partial x^2}{\partial \bar{x}} &= \frac{\bar{x}}{\bar\rho\bar{R}^2}\left( C_2^z-C_1^z + R_2\cos\theta_2-R_1\cos\theta_1\right),\\ \frac{\partial x^2}{\partial \bar{y}} &= \frac{\bar{y}}{\bar\rho\bar{R}^2}\left( C_2^z-C_1^z + R_2\cos\theta_2-R_1\cos\theta_1\right),\\ \frac{\partial x^2}{\partial \bar{z}} &= \frac{1}{2}(1-\lambda)(z^{+}_{\mathrm{P}1}-z^{-}_{\mathrm{P}1})+ \frac{1}{2}\lambda(z^{+}_{\mathrm{P}2}-z^{-}_{\mathrm{P}2}). \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
\begin{align} 0.98 R_2 &\geq R_1 + |C_1-C_2|, \label{eq:spherecontained} \end{align}
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 \(R_1 \geq 0.08 R_2\). Again, this assumption is made for accuracy purposes and might be relaxed.
Consider the line segment \(L^+_1\) 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^+_2\) 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^+_1\) and \(L^+_2\), 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{min}}\), then the line segment \(L^+_2\) twice intersects the unmapped portion of sphere 1 near the north pole, so the map is ill-defined. Similarly, if \(\alpha^+ < \theta_{2 \mathrm{min}}\), then the line segment \(L^+_2\) twice intersects the mapped portion of sphere 2 near the north pole, and again the map is poorly defined. Therefore we demand that the map parameters satisfy
where 1.1 is a safety factor.
Similarly, one can define an angle \(\alpha^-\) for the region near the south pole, and we require similar restrictions on that angle.
We also demand that either \(z^+_{\mathrm{P}1} = z^+_{\mathrm{P}2}\) or that \(z^+_{\mathrm{P}1} <= z^+_{\mathrm{P}2} -0.03 R_2\). Similarly, we demand that either \(z^-_{\mathrm{P}1} = z^-_{\mathrm{P}2}\) or \(z^-_{\mathrm{P}1} >= z^-_{\mathrm{P}2} + 0.03 R_2\). These restrictions follow expected use cases and avoid extreme distortions.
For \(z^+_{\mathrm{P}1} \neq z^+_{\mathrm{P}2}\) and \(z^-_{\mathrm{P}1} \neq z^-_{\mathrm{P}2}\), we assume the following restrictions on other parameters:
We prohibit a tiny Sphere 1 near the edge of Sphere 2 by demanding that
\begin{align} C^z_1 - R_1 &\leq C^z_2 + R_2/5,\\ C^z_1 + R_1 &\geq C^z_2 - R_2/5. \end{align}
We also demand that the polar axis of Sphere 2 intersects Sphere 1 somewhere:
\begin{align} \sqrt{(C^x_1-C^x_2)^2 + (C^y_1-C^y_2)^2} &\leq R_1. \end{align}
and we demand that Sphere 1 is not too close to the edge of Sphere 2 in the \(x\) or \(y\) directions:
\begin{align} \sqrt{(C^y_1-C^y_2)^2 + (C^y_1-C^y_2)^2} &\leq \mathrm{max}(0,0.95 R_2-R_1), \end{align}
where the max avoids problems when \(0.95 R_2-R_1\) is negative (which, if it occurs, means that the \(x\) and \(y\) centers of the two spheres are equal).
We require that the z planes in the above figures lie above/below the centers of the corresponding spheres and are not too close to the centers or edges of those spheres; specificially, we demand that
\begin{align} \label{eq:theta_1_min_res} 0.15\pi &< \theta_{1 \mathrm{min}} < 0.4\pi \\ \label{eq:theta_1_max_res} 0.6\pi &< \theta_{1 \mathrm{max}} < 0.85\pi \\ \label{eq:theta_2_min_res} 0.15\pi &< \theta_{2 \mathrm{min}} < 0.4\pi \\ \label{eq:theta_2_max_res} 0.6\pi &< \theta_{2 \mathrm{max}} < 0.85\pi . \end{align}
Here the numerical values are safety factors. These restrictions are not strictly necessary but are made for simplicity. Increasing the range will make the maps less accurate because the domain is more distorted. These parameters can be changed provided the unit tests are changed to test the appropriate parameter ranges.
If \(z^+_{\mathrm{P}1} = z^+_{\mathrm{P}2}\) or \(z^-_{\mathrm{P}1} = z^-_{\mathrm{P}2}\) we demand that \(C_1^x=C_2^x\) and \(C_1^y=C_2^y\), which simplifies the cases we need to test and agrees with our expected use cases. We also demand
\begin{align} z^+_{\mathrm{P}2} &\geq z^-_{\mathrm{P}2} + 0.18 R_2 \end{align}
This condition is necessary because for unequal z planes, \(\theta_{2 \mathrm{min}}\) and \(\theta_{2 \mathrm{max}}\) are no longer required to be on opposite sides of the equator of sphere 2 (see the paragraph below). Note that for unequal z planes \(\theta_{1 \mathrm{min}}\) and \(\theta_{1 \mathrm{max}}\) are no longer required to be on opposite sides of the equator of sphere 1, but the conditions in the paragraph below guarantee that \(z^+_{\mathrm{P}1} \geq z^-_{\mathrm{P}1}\).
Unlike the case with unequal z planes, we no longer require that the z planes in the above figures lie above/below the centers of the corresponding spheres, but we still require that the z planes are not too close to the edges of those spheres. The restrictions are the same as Eqs. ( \(\ref{eq:theta_1_min_res}\)– \(\ref{eq:theta_2_max_res}\)) except for the following changes: If \(z^+_{\mathrm{P}1} = z^+_{\mathrm{P}2}\), then we replace Eq. ( \(\ref{eq:theta_1_min_res}\)) with
\begin{align} \label{eq:equal_plus_theta_1_min_res} 0.15\pi &< \theta_{1 \mathrm{min}} < 0.59\pi, \end{align}
and furthermore, if \(z^+_{\mathrm{P}1} = z^+_{\mathrm{P}2}\) and \(\theta_{1 \mathrm{min}} > 0.4\pi\) we replace Eqs. ( \(\ref{eq:theta_1_max_res}\)– \(\ref{eq:theta_2_min_res}\)) with
\begin{align} \label{eq:equal_plus_high_theta_1_max_res} 0.7\pi &< \theta_{1 \mathrm{max}} < 0.85\pi \\ \label{eq:equal_plus_high_theta_2_min_res} 0.25\pi &< \theta_{2 \mathrm{min}} < 0.75\pi, \end{align}
but if \(z^+_{\mathrm{P}1} = z^+_{\mathrm{P}2}\) and \(\theta_{1 \mathrm{min}} \leq 0.4\pi\) we replace Eq. ( \(\ref{eq:theta_2_min_res}\)) with
\begin{align} \label{eq:equal_plus_low_theta_2_min_res} 0.15\pi &< \theta_{2 \mathrm{min}} < 0.75\pi. \end{align}
Similarly, if \(z^-_{\mathrm{P}1} = z^-_{\mathrm{P}2}\) we replace ( \(\ref{eq:theta_1_max_res}\)) with
\begin{align} \label{eq:equal_minus_theta_1_max_res} 0.41\pi &< \theta_{1 \mathrm{max}} < 0.85\pi, \end{align}
and furthermore, if \(z^-_{\mathrm{P}1} = z^-_{\mathrm{P}2}\) and \(\theta_{1 \mathrm{max}} < 0.6\pi\) we replace Eqs. ( \(\ref{eq:theta_1_min_res}\)) and ( \(\ref{eq:theta_2_max_res}\)) with
\begin{align} \label{eq:equal_minus_high_theta_1_min_res} 0.15\pi &< \theta_{1 \mathrm{min}} < 0.3\pi \\ \label{eq:equal_minus_high_theta_2_max_res} 0.25\pi &< \theta_{2 \mathrm{max}} < 0.75\pi, \end{align}
but if \(z^-_{\mathrm{P}1} = z^-_{\mathrm{P}2}\) and \(\theta_{1 \mathrm{max}} \geq 0.6\pi\) we replace Eq. ( \(\ref{eq:theta_2_max_res}\)) with
\begin{align} \label{eq:equal_minus_low_theta_2_max_res} 0.25\pi &< \theta_{2 \mathrm{max}} < 0.85\pi . \end{align}