As a first step, I'd move $P_1$ to the origin, so that the points become $P_1=(0,0,0)$, $P_2=(x_2-x_1, y_2-y_1, z_2-z_1)$, $P_3=(x_3-x_1, y_3-y_1, z_3-z_1)$. From there it's just a question of applying a rotation matrix.
A rotation matrix $\mathbf{R}$ which preserves all distances and shapes etc is
given by
$$
\mathbf{R}=\left(
\begin{array}
[c]{ccc}%
e_{1}^{x} & e_{2}^{x} & e_{3}^{x}\\
e_{1}^{y} & e_{2}^{y} & e_{3}^{y}\\
e_{1}^{z} & e_{2}^{z} & e_{3}^{z}%
\end{array}
\right)
$$
where the vectors $\left( e_{1}^{\alpha},e_{2}^{\alpha},e_{3}^{\alpha
}\right) $ for $\alpha\in\left\{ x,y,z\right\} $ are orthogonal and of unit
length, i.e.
$$
\sum_{i=1}^{3}e_{i}^{\alpha}e_{i}^{\beta}=\left\{
\begin{array}
[c]{ccc}%
1 & & \alpha=\beta\\
0 & & \alpha\neq\beta
\end{array}
\right. .%
$$
Applying $\mathbf{R}$, your rotated points $P_{i}^{\prime}=\left(
x_{i}^{\prime},y_{i}^{\prime},z_{i}^{\prime}\right) $ are given by
$$
\left(
\begin{array}
[c]{c}%
x_{i}^{\prime}\\
y_{i}^{\prime}\\
z_{i}^{\prime}%
\end{array}
\right) =\mathbf{R}\left(
\begin{array}
[c]{c}%
x_{i}\\
y_{i}\\
z_{i}%
\end{array}
\right)
$$
I suggest making the top row of $\mathbf{R}$ the unit
vector in the direction from $P_{1}$ to $P_{2}$, i.e.
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{x}\\
e_{2}^{x}\\
e_{3}^{x}%
\end{array}
\right) =\frac{1}{\sqrt{\left( x_{1}-x_{2}\right) ^{2}+\left( y_{1}%
-y_{2}\right) ^{2}+\left( z_{1}-z_{2}\right) ^{2}}}\left(
\begin{array}
[c]{c}%
x_{1}-x_{2}\\
y_{1}-y_{2}\\
z_{1}-z_{2}%
\end{array}
\right)
$$
which would mean that the line from $P_{1}$ to $P_{2}$ gets mapped into the
$x$-axis.
Working out what you could use for the second and third rows of $\mathbf{R}$ is now just a question of solving some simple linear equations, to make sure that the line from $P_{1}$ to $P_{3}$ has no z component.
To solve for $\left( e_{1}^{z},e_{2}^{z},e_{3}^{z}\right)$, the vector
$e_{i}^{z}$ must have a zero dot-product with both $e_{i}^{x}$ and $\left(
x_{3}-x_{1},y_{3}-y_{1},z_{3}-z_{1}\right) $, so the equations are
$$
\begin{align*}
e_{1}^{x}e_{1}^{z}+e_{2}^{x}e_{2}^{z}+e_{3}^{x}e_{3}^{z} & =0\\
\left( x_{3}-x_{1}\right) e_{1}^{z}+\left( y_{3}-y_{1}\right) e_{2}%
^{z}+\left( z_{3}-z_{1}\right) e_{3}^{z} & =0
\end{align*}
$$
Hence
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{z}\\
e_{2}^{z}\\
e_{3}^{z}%
\end{array}
\right) =\lambda_{z}\left(
\begin{array}
[c]{c}%
\left( z_{3}-z_{1}\right) e_{2}^{x}-\left( y_{3}-y_{1}\right) e_{3}^{x}\\
\left( x_{3}-x_{1}\right) e_{3}^{x}-\left( z_{3}-z_{1}\right) e_{1}^{x}\\
\left( y_{3}-y_{1}\right) e_{1}^{x}-\left( x_{3}-x_{1}\right) e_{2}^{x}%
\end{array}
\right)
$$
for some $\lambda_z$, which should be determined so that $e_{i}^{z}$ is a vector
of unit length. Finally $e_{i}^{y}$ can be determined as the vector which is
perpendicular to both $e_{i}^{x}$ and $e_{i}^{z}$, i.e.
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{y}\\
e_{2}^{y}\\
e_{3}^{y}%
\end{array}
\right) =\lambda_{y}\left(
\begin{array}
[c]{c}%
e_{3}^{z}e_{2}^{x}-e_{2}^{z}e_{3}^{x}\\
e_{1}^{z}e_{3}^{x}-e_{3}^{z}e_{1}^{x}\\
e_{2}^{z}e_{1}^{x}-e_{1}^{z}e_{2}^{x}%
\end{array}
\right)
$$
where again $\lambda_{y}$ is determined so that $e_{i}^{y}$ is a vector of
unit length.
The points on the infinite right circular cylinder are at distance $r$ from the axis of the cylinder. The distance between point $\vec{p}$ and line through $\vec{p}_1$ and $\vec{p}_2$ is
$$d = \frac{\left \lvert (\vec{p} - \vec{p}_1) \times (\vec{p} - \vec{p}_2) \right \rvert}{\left \lvert \vec{p}_1 - \vec{p}_2 \right \rvert}$$
so we can square that to get our equation describing an infinite right circular cylinder of radius $r$ whose axis passes through points $\vec{p}_1$ and $\vec{p}_2$:
$$r^2 = \frac{\left \lvert (\vec{p} - \vec{p}_1) \times (\vec{p} - \vec{p}_2) \right \rvert^2}{\left \lvert \vec{p}_1 - \vec{p}_2 \right \rvert^2}$$
That can be written using Cartesian coordinates as
$$r^2 = \frac{ \left ( (y-y_1)(z-z_2) - (z-z_1)(y-y_2) \right )^2 }{(x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1)^2 } + \frac{ \left ( (z-z_1)(x-x_2) - (x-x_1)(z-z_2) \right )^2 }{(x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1)^2 } + \frac{ \left ( (x-x_1)(y-y_2) - (y-y_1)(x-x_2) \right )^2 }{(x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1)^2 }$$
which we can convert to standard form by multiplying by the right-hand-side denominator, and moving all to the same side:
$$\left ( (y-y_1)(z-z_2) - (z-z_1)(y-y_2) \right )^2 + \left ( (z-z_1)(x-x_2) - (x-x_1)(z-z_2) \right )^2 + \left ( (x-x_1)(y-y_2) - (y-y_1)(x-x_2) \right )^2 - r^2 \left ( (x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1)^2 \right ) = 0$$
I have at my Wikipedia talk page the equations to locate the intersection between a line passing through origin and an arbitrary right circular cylinder; without caps, with flat end caps, or with spherical end caps. The key point is that when the ray (or line we are testing intersection for) passes through origin (so we use an unit vector $\hat{n}$, $\lvert\hat{n}\rvert=1$, to describe the line or ray), the equations become rather straightforward.
If the underlying problem is related to grinding spheres using a grinding wheel -- a cutting plane, essentially -- then a simplified coordinate system where origin is the sphere origin, and $\hat{z}$ is along the rotation axis, makes for easy calculations. Consider this illustration:
![Grinding sphere]](https://web.archive.org/web/20190515124011/http://www.nominal-animal.net/answers/tilted-sphere.svg)
The contact point scribes a circle around the axis of rotation. We can describe this circle using a single scalar, $d$, which is the signed distance from the origin to the plane of that circle, along the axis of rotation. If the circle is not on the surface of the sphere, but perhaps inside it, then we need a second scalar, $r$, to describe the radius of the circle (around the axis of the rotation).
If the radius of the sphere is $R$, then on the surface of the sphere $r = \sqrt{R^2-d^2}$, of course.
If we include the cutting or abrasion done by the wheel, then two pairs of scalars are needed: $(d_1,r_1)$ and $(d_2,r_2)$. If we assume the grinding wheel is planar, and stays at a fixed location with a perfect cut for at least one rotation of the sphere, the removed section of the sphere leaves a flat facet. Mathematically, $$r(d) = r_1 + \left ( d - d_1 \right ) \frac{r_2 - r_1}{d_2 - d_1}$$
where $d_1 \le d \le d_2$. Remember, $d$ is the distance along the axis of rotation, and $r$ the circle radius (around a point on the axis of rotation, at distance $d$ from the center of the sphere); that is why the dependence above is linear too.
If the axis of the grinding wheel always intersects the rotation axis, and you call up $z$ and right $x$ in the above diagram, then $(d,r) = (z,x)$, and the primary problem (surface of revolution caused by grinding, or shape of the "sphere" after several grinds) reduces to planar (2D) vector calculus.
If we also have the axis of rotation $\hat{n}$ in some fixed sphere coordinate system, it is easy to calculate the "ribbons" each grind cuts into the sphere, either parametrically (using say $0 \le u \le \pi$ along the circular ribbon, and $0 \le v \le 1$ across the ribbon -- useful for visualization) or algebraically.
![Grinding sphere]](http://www.nominal-animal.net/answers/tilted-sphere.svg)
Best Answer
We can find the distance $d$ from point $\vec{p}$ to the infinite right circular cylinder whose axis passes through points $\vec{p}_1$ and $\vec{p}_2$ using point-line distance (to the axis): $$d = \frac{ \left\lVert \left ( \vec{p}_2 - \vec{p}_1 \right ) \times \left ( \vec{p}_1 - \vec{p} \right ) \right\rVert}{\left\lVert \vec{p}_2 - \vec{p}_1 \right\rVert}$$
We can use $t$ to denote the position along the axis (so that the line between that point and point $\vec{p}$ is perpendicular to the axis), with $t=0$ at $\vec{p}_1$, and $t=1$ at $\vec{p}_2$: $$t = \frac{ \left( \vec{p} - \vec{p}_1 \right) \cdot \left( \vec{p}_2 - \vec{p}_1 \right) }{\left\lVert \vec{p}_2 - \vec{p}_1 \right\rVert^2}$$
Point $\vec{p}$ is within the right circular cylinder, if and only if $$\begin{cases} t \ge 0 \\ t \le 1 \\ d \le r \end{cases}$$
Because the right circular cylinder is axially symmetric, there is no "zero $\theta$ angle" we can refer to. We must establish one ourselves.
Let's say we have an unit vector $\hat{c}$ ($\left\lVert\hat{c}\right\rVert^2 = \hat{c}\cdot\hat{c} = 1$) that establishes that zero angle for us. The vector must be perpendicular to the axis, i.e. $$\hat{c} \cdot \left ( \vec{p}_2 - \vec{p}_1 \right ) = 0$$ If you have some generic zero angle vector $\vec{v}$ which is not parallel to the axis, you can calculate $\hat{c}$ in two steps, using $$\vec{v}' = \vec{v} - \vec{v} \cdot \left( \vec{p}_2 - \vec{p}_1 \right) \frac{\vec{p}_2 - \vec{p}_1}{\left\lVert \vec{p}_2 - \vec{p}_1 \right\rVert^2}$$ Note that $\vec{v}'$ is the original vector $\vec{v}$ with the axis-parallel-part subtracted from it; thus, it is perpendicular to the axis. We only need to normalize it to unit length to get $\hat{c}$: $$\hat{c} = \frac{\vec{v}'}{\left\lVert\vec{v}'\right\rVert}$$
If we already know $d$ and $t$, we can easily construct an unit vector $\hat{q}$, from the axis towards point $\vec{p}$: $$\hat{q} = \frac{\vec{p} + t \left ( \vec{p}_1 + \vec{p}_2 \right ) - \vec{p}_1}{\left\lVert \vec{p} + t \left ( \vec{p}_1 + \vec{p}_2 \right ) - \vec{p}_1 \right\rVert}$$ This unit vector is perpendicular to the axis (because that's how $t$ was defined), so the angle $\theta$ can now be easily calculated: $$\begin{cases} \cos\theta = \hat{q} \cdot \hat{c} \\ \sin\theta = \left\lVert \hat{q} \times \hat{c} \right\rVert \end{cases}$$
When $\vec{p}_1$, $\vec{p}_2$, $\hat{c}$, $t$, $d$, $\sin\theta$, and $\cos\theta$ are known, we can calculate point $\vec{p}$.
First, we need a helper unit vector $\hat{u}$ which is perpendicular to both the axis ($\hat{u}\cdot\left(\vec{p}_2-\vec{p}_1\right)=0$) and the zero angle unit vector $\hat{c}$ ($\hat{u}\cdot\hat{c}=0$): $$\hat{u} = \frac{ \hat{c} \times \left ( \vec{p}_2 - \vec{p}_1 \right ) }{\left\lVert \hat{c} \times \left ( \vec{p}_2 - \vec{p}_1 \right ) \right\rVert }$$
Then, $$\vec{p} = (1 - t)\,\vec{p}_2 + t\,\vec{p}_1 + \hat{c}\,d \cos\theta + \hat{u}\,d \sin\theta$$
Years ago, I wrote some related equations to my Wikipedia talk page, but never found the effort to try and push it anywhere more approppriate to be worthwhile. You might find it interesting or useful, though.