Intersection of sphere and cone with axis along arbitrary direction I need to find the coordinates of the circle formed by intersecting a sphere of radius $R$ centred at the origin with a cone of opening half-angle $\alpha$ that has its vertex at the origin and axis pointing along an arbitrary direction given by the unit vector $\mathbf{\hat{n}} = (\sin \theta_C \,\cos \phi_C, \;\sin \theta_C \,\sin \phi_C, \;\cos \theta_C)$.
I've taken the cone to be described as the set of vectors that form an angle $\alpha$ with the axis of the cone (I'm really only interested in the half of the cone that points along the direction of the unit vector, not anti-parallel to it). So I took
$$\mathbf{x}\cdot\mathbf{\hat{n}} = \lvert\mathbf{x}\rvert \cos \alpha$$
for all vectors $\mathbf{x}$ that lie on the cone. In terms of Cartesian coordinates, this becomes
$$prx + psy + qz = \sqrt{x^2+y^2+z^2} \cos\alpha$$
where $p = \sin \theta_C,\;q=\cos \theta_C,\;r=\cos \phi_C,\;s=\sin \phi_C$ to keep it clean.
For the sphere, I've got $$x^2+y^2+z^2=R^2$$
Intersecting the cone and the sphere by expressing $z=\pm\sqrt{R^2-x^2-y^2}$ and substituting, I get
$$\pm q \sqrt{R^2-x^2-y^2}=R \cos\alpha - prx -psy$$ and upon squaring
$$x^2(p^2r^2+q^2) + y^2(p^2s^2+q^2) + 2xyp^2rs - 2R\cos\alpha\,p(xr+ys) = R^2 (q^2-cos^2\alpha)$$
From this, I solved for y in terms of x and then substituted for z to get two branches of what should be a circle. (The reason I did not do this in spherical coordinates is that I need to find the "start" and "end" points of the circle in the independent variable and with $\sin$ and $\cos$ corresponding to more than one possible angle, this might be ambiguous.) I think so far this should be correct.
My question is: Since the intersection in this configuration should be a circle of radius $r = R \sin\alpha$, I should be able to express the above equation in the form $$(x')^2 + (y')^2 = R^2 \sin^2 \alpha$$ for some coordinates $x',\,y'$ in the plane of the circle, right? How should I go about finding the necessary coordinate transformation to get the earlier (messy) equation in this form?
(Context: As a next step, I would like to integrate a function of spherical coordinates along this circle, so I am going to need some kind of parameterisation of the circle. I could take x as my parameter and integrate along the two branches of the circle with $y(x)$ and $z(x)$, but I'm not sure if this would be easy.)
 A: The result you're looking for is a tilted 3D circle.  If we concentrate on the means of representing this Circle3D numerically, the answer comes pretty easily.
Desired result:   Circle3D representation (C-O-R):
Center  (point vector)
Orientation (direction vector)
Radius (scalar)
Problem inputs:
Sphere SPH =  [ (0 , 0, 0 ] , r ] 
Cone CONE = [ ( 0, 0, 0 ), d_axis, half_angle ]
Solution:
The circle Center is on d_axis, and at distance rcos(half_angle) from the origin, so a scalar multiply locates Center:
Center <--  CONE.d_axis.scalarMult( r * cos(CONE.half-angle) )
The result circle will have orientation defined by CONE's axis direction:
Orientation <-- CONE.d_axis
All that's left to solve is the circle radius.  We know the half-arc length (along the sphere surface) from the cone axis out to the circle edge:
half_arc_length = r * CONE.half_angle     (half_angle in radians)
The circle radius is the half-chord length (straight-line distance between two endpoints of the half-arc:
Radius = r * sin (CONE.half_angle)
A: Since you want to end up with a parameterization of the sphere-cone intersection, instead of trying to manipulate the messy equation further, you can start with the parameterized circle $(R\sin\alpha\cos t,R\sin\alpha\sin t)$ and look for a transformation that maps this onto the cone-sphere intersection. This transformation will consist of a rotation and translation.$$ Taking advantage of the fact that the columns of a transformation matrix are the images of the basis vectors, this problem reduces to finding the a suitable orthonormal basis for $\mathbb R^3$.  
The translation part is easy: the center of the circle gets shifted from the origin to $\mathbf c=(R\cos\alpha)\mathbf n$, which is the intersection of the cone’s axis with the circle’s plane.
For the rotation part, we want to map the $z$-axis to the normal of the intersection circle’s plane, which is $\mathbf n$. If $\mathbf n$ is parallel to the $z$-axis, no rotation is necessary—the translation is sufficient by itself, we can take the standard basis for our orthonormal basis and we’re done. Otherwise, we need to find an orthonormal pair of unit vectors parallel to the intersection circle’s plane. This can be done via cross products. If $\mathbf n=(n_x,n_y,n_z)^T$, then let $$\mathbf u={\mathbf n\times(0,0,1)^T\over\|\mathbf n\times(0,0,1)^T\|} = {(n_y,-n_x,0)^T\over\sqrt{n_x^2+n_y^2}}$$ and $\mathbf v=\mathbf n\times\mathbf u$. Be careful with the order of the vectors in these products so that you don’t change the orientation of the circle.  
With these two vectors in hand, we can construct the parameterization we seek: $$\mathbf c+R\sin\alpha\cos t\,\mathbf u+R\sin\alpha\sin t\,\mathbf v.$$ (You can derive this from the original parameterization by constructing the transformation matrix and applying it to the homogeneous vector $[R\sin\alpha\cos t,R\sin\alpha\sin t,0,1]^T$.) Any pair of orthonormal vectors that are orthogonal to $\mathbf n$ will do for this parameterization: Choosing a different vector for $\mathbf u$ moves the “start” point of the circle.
A: This may not be the answer you want, but I suspect it's the answer that you should want, so I'm gonna write it. 
Step 1: Find a rotation that takes the north pole to your vector $n$. (There are infinitely many, but pick one). Represent this by multiplication by a matrix $M$. 
Step 2: Parameterize the circle with cone-angle $\alpha$ around the north pole. If the north pole is $(0,0,1)$ (I like the  "z is up" convention rather thna the "y is up" convention) then this is something like 
$$
u(t) = \pmatrix{
R\sin \alpha \sin t\\
R \sin \alpha \cos t\\
R \cos \alpha}
$$
Step 3: Rotate $u$ into the position you want:
$$
v(t) = M \cdot u(t).
$$
For step 1, you can form the vector $w' = (0,0,1) \times n$, and then normalize it (and call the result $w$). Then compute $\theta = \cos^{-1} (n \cdot (0,0,1))$, which is what I think you've called $\theta_C$. Now use Rodrigues' formula to build a rotation by angle $-\theta$ around $w$. (It may have to be $\theta$...try it once to find out!)
This part-1 trick fails if $n$ is the north pole...but in that case, you don't need the rotation part, and can choose it to be the identity matrix.  
