Solution (source code)

= Solution

For $v\in T_pS$ for which the <geodesic> exists up to time one, the <exponential map> is $\exp_p(v)=\gamma_v(1)$, where $\gamma_v(0)=p$ and $\dot\gamma_v(0)=v$. The initial velocity determines the <geodesic> uniquely. Since $d(\exp_p)_0$ is the identity, this map is a local <diffeomorphism> near zero.

Choose an orthonormal <basis> of $T_pS$ and set $e(\theta)=(\cos\theta,\sin\theta)$. The <geodesic polar coordinates> are
$$
\Phi(r,\theta)=\exp_p(re(\theta))
$$
in a normal neighbourhood, away from its centre. For sufficiently small positive $r$, the <geodesic circle> is
$$
\boxed{S_r(p)=\{\exp_p(re(\theta)):0\le\theta<2\pi\}
=\{x:d(p,x)=r\}}
$$
within that neighbourhood. Outside the injectivity range the distance sphere need not have this smooth coordinate description.