Solution (source code)

= Solution

In <Cartesian coordinates>, $\nabla\Phi=(\Phi_x,\Phi_y,\Phi_z)$. For <orthogonal coordinates>, let $\mathbf r(q_1,q_2,q_3)$ be the position, $h_i=|\partial\mathbf r/\partial q_i|$ the <scale factors of orthogonal coordinates>, and $\mathbf e_i=h_i^{-1}\partial\mathbf r/\partial q_i$ the corresponding <unit vectors>. Since $d\mathbf r=\sum_i h_i\mathbf e_i\,dq_i$ and $d\Phi=\nabla\Phi\cdot d\mathbf r=\sum_i\Phi_{q_i}dq_i$, comparison of independent increments gives the <gradient in orthogonal coordinates>:
$$
\boxed{\nabla\Phi=\sum_{i=1}^3\frac{\mathbf e_i}{h_i}\frac{\partial\Phi}{\partial q_i}.}
$$
For the <hyperbolic parametrization of oblate spheroidal coordinates>, the three coordinate tangent vectors are
$$
\mathbf r_R=(\sinh R\cos\theta\cos\phi,\sinh R\cos\theta\sin\phi,\cosh R\sin\theta),
$$
$$
\mathbf r_\theta=(-\cosh R\sin\theta\cos\phi,-\cosh R\sin\theta\sin\phi,\sinh R\cos\theta),\qquad\mathbf r_\phi=(-\cosh R\cos\theta\sin\phi,\cosh R\cos\theta\cos\phi,0).
$$
Their pairwise <dot products> are zero. Their squared <Euclidean norms> are $\sinh^2R+\sin^2\theta$, $\sinh^2R+\sin^2\theta$, and $\cosh^2R\cos^2\theta$, respectively. Thus the coordinate surfaces meet orthogonally wherever this is a regular coordinate chart, and
$$
\boxed{h_R=h_\theta=\sqrt{\sinh^2R+\sin^2\theta},\qquad h_\phi=\cosh R\cos\theta.}
$$
The last expression uses the usual latitude range $|\theta|<\pi/2$; without that range a scale factor is its absolute value. Axial or ring degeneracies are excluded from a regular chart.