Circular arc under Lambert Azimuthal Equal Area projection I wish to produce a formula to project an arc of a great circle on a unit sphere.  The arc is simply defined by 2 points (x, y, z), and goes the shortest distance between them.  The LAEA projection takes (x, y, z) to $k=\sqrt{2/(1-z)}$, $u=x/k$, $v=y/k$.  I would like to use something like a cubic bézier curve to draw the projected arc, to good approximation.  
My 1st idea was to construct the bézier curve that matches the circular arc, which is fairly easy to do.  Then I attempted to project the control points of that curve onto the (u, v) plane, and draw that.  It comes close, with the endpoints obviously matching, but longer arcs deviate significantly in the middle.  I'm not sure how I need to adjust the middle 2 control points.
From the arc on the sphere, I can produce any number of points on the projection plane.  From the bézier control points on the sphere, I have vectors from the end points towards the middle points.  I'm not entirely sure how to convert those 3D vectors into useful 2D vectors on the projection plane (since they aren't on the sphere, and the LAEA projection is really only meaningful from the sphere), or how I could easily use those vectors to find the middle control points on the plane.
Any suggestions would be helpful.
[Edit: By the way, it's quite alright if the arc is limited to under 90 degrees, but solutions should be valid for endpoints anywhere, except possibly (0,0,1) (division by 0 on $k$).]

If the endpoints of your arc are $p_{1} = (x_{1}, y_{1}, z_{1})$ and $p_{2} = (x_{2}, y_{2}, z_{2})$, then the angle between these points (measured from the center of the sphere) is
$$
\theta = \arccos (p_{1} \cdot p_{2})
  = \arccos(x_{1}x_{2} + y_{1}y_{2} + z_{1}z_{2}).
$$
The unit vector
$$
u = \frac{p_{2} - (p_{1} \cdot p_{2}) p_{1}}{\|p_{2} - (p_{1} \cdot p_{2}) p_{1}\|}
$$
(obtained from the Gram-Schmidt algorithm, if it's of interest) lies in the plane containing $p_{1}$ and $p_{2}$, and is orthogonal to $p_{1}$. The (short) great circle arc from $p_{1}$ to $p_{2}$ is parametrized by
$$
(x, y, z) = (\cos t) p_{1} + (\sin t)u,\qquad 0 \leq t \leq \theta.
\tag{1}
$$
Using (1), you can generate as many points exactly on the arc as you like and map them to the plane, then use them as control points.
Geometrically, this scheme amounts to looking at the plane through the origin, $p_{1}$, and $p_{2}$; introducing Cartesian coordinates so that $p_{1}$ corresponds to $e_{1} = (1, 0)$ and $p_{2}$ corresponds to $(\cos\theta, \sin\theta)$ with $0 < \theta < \pi$; calculating the spatial vector $u$ corresponding to $e_{2} = (0, 1)$; and parametrizing the arc as $(\cos t) e_{1} + (\sin t)e_{2}$ for $0 \leq t \leq \theta$.
The computation of $u$ presumes the points $p_{1}$ and $p_{2}$ are distinct and not antipodal, and becomes numerically unstable if $p_{1}$ and either of $\pm p_{2}$ are close to each other.
