cartesian_to_lat_lon(x, y, z)

Convert 3D cartesian coordinates (x, y, z) on the sphere to latitude-longitude. Returns a tuple (latitude, longitude) in degrees.

The equatorial plane falls at $z = 0$, latitude is the angle measured from the equatorial plane, and longitude is measured anti-clockwise (eastward) from $x$-axis ($y = 0$) about the $z$-axis.


Find latitude-longitude of the North Pole

julia> using CubedSphere

julia> x, y, z = (0, 0, 6.4e6); # cartesian coordinates of North Pole [in meters]

julia> cartesian_to_lat_lon(x, y, z)
(90.0, 0.0)

Let's confirm that for few points on the unit sphere we get the answers we expect.

julia> cartesian_to_lat_lon(√2/4, -√2/4, √3/2)
(59.99999999999999, -45.0)

julia> cartesian_to_lat_lon(-√6/4, √2/4, -√2/2)
(-45.00000000000001, 150.0)
conformal_cubed_sphere_inverse_mapping(X, Y, Z; Z_map=Z_Rancic)

Inverse mapping for conformal cube sphere for quadrant of North-pole face in which X and Y are both positive. All other mappings to other cube face coordinates can be recovered from rotations of this map. There a 3 other quadrants for the north-pole face and five other faces for a total of twenty-four quadrants. Because of symmetry only the reverse for a single quadrant is needed. Because of branch cuts and the complex transform the inverse mappings are multi-valued in general, using a single quadrant case allows a simple set of rules to be applied.

The mapping is valid for the cube face quadrant defined by $0 < x < 1$ and $0 < y < 1$, where a full cube face has extent $-1 < x < 1$ and $-1 < y < 1$. The quadrant for the mapping is from a cube face that has "north-pole" at its center $(x=0, y=0)$. i.e., has X, Y, Z = (0, 0, 1) at its center. The valid ranges of X and Y for this mapping and convention are a quadrant defined be geodesics that connect the points A, B, C and D, on the shell of a sphere of radius $R$ with X, Y coordinates as follows

A = (0, 0)
B = (√2, 0)
C = (√3/3, √3/3)
D = (0, √2)
conformal_cubed_sphere_mapping(x, y; W_map=W_Rancic)

Conformal mapping from a face of a cube onto the equivalent sector of a sphere with unit radius.

Map the north-pole face of a cube with coordinates $(x, y)$ onto the equivalent sector of the sphere with coordinates $(X, Y, Z)$.

The cube's face oriented normal to $z$-axis and its coordinates must lie within the range $-1 ≤ x ≤ 1$, $-1 ≤ y ≤ 1$ with its center at $(x, y) = (0, 0)$. The coordinates $X, Y$ increase in the same direction as $x, y$.

The numerical conformal mapping used here is described by Rancic-etal-1996.

This is a Julia translation of MATLAB code from MITgcm that is based on Fortran 77 code from Jim Purser & Misha Rančić.


The center of the cube's face $(x, y) = (0, 0)$ is mapped onto $(X, Y, Z) = (0, 0, 1)$

julia> using CubedSphere

julia> conformal_cubed_sphere_mapping(0, 0)
(0, 0, 1.0)

and the edge of the cube's face at $(x, y) = (1, 1)$ is mapped onto $(X, Y, Z) = (√3/3, √3/3, √3/3)$

julia> using CubedSphere

julia> conformal_cubed_sphere_mapping(1, 1)
(0.5773502691896256, 0.5773502691896256, 0.5773502691896257)


find_N(r; decimals=15)

Return the required number of points we need to consider around the circle of radius r to compute the conformal map series coefficients up to decimals points. The number of points is computed based on the estimate of eq. (B9) in the paper by Rancic-etal-1996. That is, is the smallest integer $N$ (which is chosen to be a power of 2 so that the FFTs are efficient) for which

\[N - \frac{7}{12} \frac{\mathrm{log}_{10}(N)}{\mathrm{log}_{10}(r)} - \frac{r + \mathrm{log}_{10}(A₁ / C)}{-4 \mathrm{log}_{10}(r)} > 0\]

where $r$ is the number of decimals we are aiming for and

\[C = \frac{\sqrt{3} \Gamma(1/3) A₁^{1/3}}{256^{1/3} π}\]

with $A₁$ an estimate of the $Z^1$ Taylor series coefficient of $W(Z)$.

For $A₁ ≈ 1.4771$ we get $C ≈ 0.265$.


find_taylor_coefficients(r = 1 - 1e-7;
                         Niterations = 30,
                         maximum_coefficients = 256,
                         Nevaluations = find_N(r; decimals=15))

Return the Taylor coefficients for the conformal map $Z \to W$ and also of the inverse map, $W \to Z$, where $Z = z^4$ and $W = w^3$. In particular, it returns the coefficients $A_k$ of the Taylor series

\[W(Z) = \sum_{k=1}^\infty A_k Z^k\]

and also coefficients $B_k$ the inverse Taylor series

\[Z(W) = \sum_{k=1}^\infty B_k Z^k\]

The algorithm to obtain the coefficients follows the procedure described in the Appendix of the paper by Rancic-etal-1996


  • r (positional): the radius about the center and the edge of the cube used in the algorithm described by Rancic-etal-1996. r must be less than 1; default: 1 - 10$^{-7}$.

  • maximum_coefficients (keyword): the truncation for the Taylor series; default: 256.

  • Niterations (keyword): the number of update iterations we perform on the Taylor coefficients $A_k$; default: 30.

  • Nevaluations (keyword): the number of function evaluations in over the circle of radius r; default find_N(r; decimals=15); see find_N.


julia> using CubedSphere

julia> using CubedSphere: find_taylor_coefficients

julia> A, B = find_taylor_coefficients(1 - 1e-4);
[ Info: Computing the first 256 coefficients of the Taylor serieses
[ Info: using 32768 function evaluations on a circle with radius 0.9999.
100.0%┣████████████████████████████████████████████┫ 30/30 [00:02<00:00, 12it/s]

julia> A[1:10]
10-element Vector{Float64}:
Reproducing coefficient table by Rančić et al. (1996)

To reproduce the coefficients tabulated by Rancic-etal-1996 use the default values, i.e., $r = 1 - 10^{-7}$.


cn(z::Complex, m::Real)

Compute the Jacobi elliptic function cn(z | m) following Abramowitz & Stegun (1964), Eq. 16.21.3.

dn(z::Complex, m::Real)

Compute the Jacobi elliptic function dn(z | m) following Abramowitz & Stegun (1964), Eq. 16.21.4.

sn(z::Complex, m::Real)

Compute the Jacobi elliptic function sn(z | m) following Abramowitz & Stegun (1964), Eq. 16.21.2.