The equation of a circle on sphere? One sphere with radius $R$, named big sphere, two points on it: $a (\text{(longitude)}_a, \text{(latitue)}_a), b(\text{(longitude)}_b, \text{(latitude)}_b)$, $\text{distance}\space (a,b)=r$, $a$ is center & $r$ is radius. There is another sphere, named little sphere. Now what is the  equation of circle of the intersection of two spheres? 
Parametric equation of the circle: $\text{longitude}=f(xx), \text{latitude}=g(xx)$
 A: This (calculating the equation of the circle intersection between two spheres) will probably help.
A: I assume $\text{distance}\ r$ is an Euclidean distance in 3D, not the length of the great circle's arc on the big sphere. Of course $0 \le r \le 2R$, with each 'equal' case causing a circle to degenerate to a single point.
Let's define a Cartesian coordinates system, so that:


*

*the origin is in the big sphere's center,

*XY plane is a plane of equator: $\varphi = 0$,

*positive latitudes $\varphi$ correspond to positive $z$,

*zero longitude $\lambda = 0$ corresponds to XZ plane $y=0$,

*longitude is measured from positive X semiaxis towards positive Y semiaxis.


Then:
$$
  \begin{cases}
    x = R \cos \varphi \cos \lambda \\
    y = R \cos \varphi \sin \lambda \\
    z = R \sin \varphi
  \end{cases}
  \quad \text{and}\quad
  \begin{cases}
    \varphi = \arcsin \frac zR \\
    \lambda = \arctan \frac yx
  \end{cases}
$$
Of course one needs to check for singularity when calculating $\lambda$ for $x=0$ and adjust the result by $+180^\circ$ when $x<0$.
Let $\rho$ = an angle between the big sphere's radii ending at $a$ and at $b$. Then $\rho = 2\arcsin \frac r{2R}$. 
Suppose $a = a_1$ = 'the North Pole': $\varphi_{a1} = 90^\circ, x=y=0, z=R$. Then the circle sought $C$ is a circle of latitude at $\varphi = 90^\circ - \rho$ (with positive latitudes corresponding to the 'northern' hemisphere).
It can be easily parameterized as
$$\begin{cases}\varphi_1 = 90^\circ - \rho \\ \lambda_1 = t\end{cases} \ \ \text{for}\ \ 0<t<360^\circ$$
and then transformed into XYZ coordinates:
$$
\begin{cases}
    x_1 = R \cos \varphi_1 \cos t \\
    y_1 = R \cos \varphi_1 \sin t \\
    z_1 = R \sin \varphi_1
\end{cases}
$$
Now suppose we shift $a$ along the Greenwich meridian to the south, to some point $a_2$ at latitude $\varphi_a$. That means the whole circle $C$ undergoes a rotation in 3D about the Y axis by an angle $\delta_\varphi = 90^\circ-\varphi_a$, measured from positive Z semiaxis towards positive X semiaxis. Each point of circle $C$ gets its coordinates tranformed like this:
$$\begin{cases}
x_2 = x_1 \cos \delta_\varphi + z_1 \sin \delta_\varphi \\
y_2 = y_1 \\
z_2 = z_1 \cos \delta_\varphi - x_1 \sin \delta_\varphi \\
\end{cases}$$
Now we need to shift $C$ along the circles of latitude so its 'center' gets into a desired position $a_3$ at longitude $\lambda_a$. This is a 3D rotation about the Z axis:
$$
\begin{cases}
x_3 = x_2 \cos \lambda_a - y_2 \sin\lambda_a \\
y_3 = y_2 \cos \lambda_a + x_2 \sin \lambda_a \\
z_3 = z_2
\end{cases}
$$
Now you only need to compose the transformations described above and possibly add another $(x,y,z)\to(\varphi,\lambda)$ transformation to get a parametric equation of the circle.
A: The circle belongs to a plane that can be found by subtracting the equations of both spheres:
$$(x^2+y^2+z^2-R^2)-((x-x_c)^2+(y-y_c)^2+(z-z_c)^2-r^2)=ax+by+cz+d=0.$$
Now plug the spherical coordinates to yield
$$aR\cos\theta\sin\phi+bR\sin\theta\sin\phi+cR\cos\phi+d=0.$$
You can express $\theta$ as a function of $\phi$ or conversely, by solving a linear trigonometric equation of the form
$$\alpha\cos t+\beta\sin t=\gamma.$$

WLOG (normalize if necessary), $\alpha=\cos u$ and $\beta=\sin u$, then $\cos(t-u)=\gamma$ gives $t$.
