Parametric equations for a great circle on the unit sphere I want to define a great circle on the unit sphere, given two points on its surface. 
The great circle is the intersection of a plane and a sphere, such that the plane touches the sphere's origin and the two points given (call them $p_1$ and $p_2$). Coordinates are Cartesian. I'm using parametrized equations in order to plug them into gnuplot.
With parameters $u_1$, $v_1$ and $u_2$, $v_2$ the two surfaces are:
$$x=cos(u_1)cos(v_1), y=cos(u_1)sin(v_1), z=sin(u_1)$$
$$x=u_2, y=v_2, z=\sqrt{{c_1^{-2}}-1}(u_2sin(c_2)+v_2cos(c_2)) $$
I'm not sure yet how $c_1$ and $c_2$ are supposed to be derived from the points' coordinates, but OK.
Solving this system of equations for the parameters I got the following:
$$u_1=arccos(\frac{u_2}{cos(v_1)})$$
$$v_1=arccos(\frac{u_2}{cos(u_1)})$$
$$u_2=cos(u_1)cos(v_1)$$
$$v_2=cos(u_1)sin(v_1)$$
Here I'm stuck; I don't know what I'm looking for. $c_1$ and $c_2$ disappeared because of the relationships I plugged in above. How can I turn these equations into the parameters that describe the great circle?
 A: If $\mathbf u$ and $\mathbf v$ are a pair of orthogonal unit vectors, then $\mathbf r(t) = \mathbf u\cos t+\mathbf v\sin t$ is a parameterization of a unit circle that lies in the plane spanned by these vectors. By definition, $p_1$ and $p_2$ are already unit vectors and if they’re not antipodes, span a unique plane. It’s convenient to take $\mathbf u = p_1$ so that $\mathbf r(0)=p_1$ and then use cross products to generate $\mathbf v$: $(p_1\times p_2)\times p_2$ lies in the span of $p_1$ and $p_2$ and is orthogonal to $p_1$. Normalize this vector and set $\mathbf v$ to the result. For this choice of $\mathbf v$, the circle will be traced toward $p_2$ along the shorter of the two arcs as $t$ increases. On the other hand, if $p_1\times p_2=0$, then the two points are either antipodes or coincide and you can take any unit vector orthogonal to $p_1$ as $\mathbf v$—the plane of the great circle is not uniquely determined.  
You could also start with the unit circle $\mathbf r(t) = (\cos t,\sin t)$ in the $x$-$y$ plane and tilt it into place, but computing the required rotation or reflection takes at least as much work as the above.
A: Your equation for the sphere is OK, but your equation for the plane excludes e.g. the plane $\{(x,y,z)|y=0\}$. What you want is the plane
$$
P = \{u_2 p_1 + v_2 p_2\} \le \mathbb R^3,
$$
where $u_2, v_2$ are the parameters. That is,
$$
x = u_2 p_1^1 + v_2 p_2^1, ~~~ y = u_2 p_1^2 + v_2 p_2^2, ~~~ z = u_2 p_1^3 + v_2 p_2^3,
$$
where the upper indices don't stand for exponents, but the components of $p_1, p_2$.
Of this plane, you can then find an orthonormal basis $w_1, w_2$ by applying the Gram‒Schmidt procedure, and then the parametric equation is
$$
u_3 \mapsto \cos(u_3) w_1 + \sin(u_3) w_2.
$$
