Modeling a line in 3D space using length and rotations rather than endpoints I'm using a cad software which can only create cylinders at the origin based on diameter d and height h. For simplicity sake, in this situation we can think of the cylinder as just a line so we can treat it as a vector. The orientation can only be changed by rotating it using rotate([x,y,z]) which rotates a certain amount of degrees around each axis and can be translated using translate([x,y,z]). I'm trying to write a function that can take two 3D points p1 and p2 and automatically create a cylinder (line in our case) using d and h then rotate and translate it so that the endpoints rest on p1 and p2.
It is easy enough for me to find the h of the line using the Pythagorean theorem, but the rotations are giving me a hard time. I used a bit of vector math I learned in a class a while back to find the angles from the line to each axis using unit vectors:
$p1 = (x_1, y_1, z_1) $
$p2 = (x_2, y_2, z_2)$
$\vec{F} = ((x_2-x_1) \hat{i} +(y_2-y_1)\hat{j} + (z_2-z_1)\hat{k})$
$\vec{u_F} = \frac{\vec{F}}{||\vec{F}||} = (\cos{\alpha} \hat{i} + \cos{\beta} \hat{j} + \cos{ \gamma} \hat{k})$
$\alpha = \arccos{\frac{x_2-x_1}{||\vec{F}||}}, \beta = \arccos{\frac{y_2-y_1}{||\vec{F}||}}, \gamma = \arccos{\frac{z_2-z_1}{||\vec{F}||}}$
This information is not the whole answer unfortunately. The angle from the " "-axis is not the same as rotated around the " "-axis. In addition, the order of the rotation changes the final result. Where do I go from here?
 A: The problem can be rephrased as follows: Rotate the axis of the cylinder, represented by a vector (say $\mathbf{a}$) into another vector, say $\mathbf{b}$. It's easy to see that in principle one could do this by rotating on the plane formed by the two vectors by an angle given by:
$$\cos\theta=\frac{\mathbf{a}\cdot\mathbf{b}}{||\mathbf{a}||~||\mathbf{b}||}$$
The axis of rotation is represented by the unit vector perpendicular to the plane which is the vector product of the two. It can be shown that:
$$\hat{\mathbf{n}}=\frac{\mathbf{a}\times\mathbf{b}}{||\mathbf{a}\times\mathbf{b||}}$$
Armed with these expressions one can calculate the rotation matrix given here. However, since the program uses only rotations along given axes, it's necessary to find angles such that the same job can be done by three subsequent rotations around the axes. Here I will assume that the x, y, z rotations are done in the order indicated.
Using the expressions for individual rotations around a particular axis, the net rotation matrix we get is:
$$R_z(\gamma)R_y(\beta)R_x(\alpha)=
\begin{pmatrix}
\cos\gamma\cos\beta &-\cos\gamma\sin\beta\sin\alpha-\sin\gamma\cos\alpha&-\cos\gamma\sin\beta\cos\alpha+\sin\gamma\sin\alpha\\
\sin\gamma\cos\beta&-\sin\alpha\sin\beta\sin\gamma+\cos\alpha\cos\gamma&-\cos\alpha\sin\beta\sin\gamma-\sin\alpha\cos\gamma\\
\sin\beta&\sin\alpha\cos\beta&\cos\alpha\cos\beta
\end{pmatrix}$$
which when equated with the rotation matrix around the axis yields:
$$\sin\beta=n_zn_x(1-\cos\theta)-n_y\sin\theta\\\tan\gamma=\frac{n_yn_x(1-\cos\theta)+n_z\sin\theta}{\cos\theta+n_x^2(1-\cos\theta)}\\
\tan\alpha=\frac{n_yn_z(1-\cos\theta)+n_x\sin\theta}{\cos\theta+n_z^2(1-\cos\theta)}$$ 
Now one can readily solve for the three required angles, making sure that all possible solutions in the interval $[0,2\pi)$ are checked to match the known rotation matrix, which is easy with a computer available.  
In the case of the problem stated above, vector $\mathbf{a}$ should be matched to the axis of the cylinder when created, and the vector $\mathbf{b}=\mathbf{p_1}-\mathbf{p_2}$.
A: Because OP is using OpenSCAD, the claims of the question aren't entirely correct: orientation can be changed not only by rotations about axes, but by other transformations as well. 
I'm going to assume that the cylinder has been rotated/translated/whatever by the OP, who can surely do this, so that one end is at the origin, and the other end is at $(0,0,1)$. (Maybe OpenSCAD aligns it along $(0,1,0)$, but I don't think that adapting to this difference will present any challenges). I'm not going to write correct openSCAD code because I haven't used it in about a year, but I'll write something that should get you going in the right direction. 
Here goes. 
There's a rotate(v, a) function in openSCAD that rotates around a vector $v$ by angle $a$. We'll use that. 
$$
\newcommand{zv}{\mathbf{\hat{z}}}
\newcommand{yv}{\mathbf{\hat{y}}}
\newcommand{sv}{\mathbf{{s}}}
$$
Let's let zvec (in code) or $\zv$ (in math) mean $[0,0,1]$, the initial axis of the cylinder. Then 


*

*Compute $v = Q - P$, the vector from the point $P$ to the point $Q$. If $v$ is zero, your input was bad (a length-zero cylinder makes no sense), and you'll have to handle this exceptional case. 

*Compute $u = \frac{v}{ \| v\|}$. In openSCAD, you do something like v = (1/norm(v))*v, I believe. If $u$ turns out to be exactly $-\zv$, then this is another special case; because the cylinder is already aligned with the $z$-axis, all you need to do is rotate it about the $x$-axis (for instance) by 180 degrees to get it where we want; you can then skip to step 5.)

*Compute $w = \frac12 (\zv + u)$, or in code, w = (0.5) * (u + zvec), and then $t = \frac{w}{\|w\|}$, or t = (1/norm(w)) * w. At this point, $t$ is a unit vector in the plane of $v$ and $\zv$, and is exactly halfway between those two vectors. 

*Rotate about $t$ by 180 degrees to take the cylinder, which used to point along the $\zv$ direction, and make it point along the $v$ direction instead. 

*Scale everything by $\| v \| $, i.e., norm(v), which is the distance from $P$ to $Q$. 

*Translate by $P$ to move the origin (one end of the cylinder) to $P$, and the other end (which is at the tip of $v = Q - P$) to the point $Q$. 

All of that stuff works mathematically, but I should caution you that the special case in step 2 is problematic: testing equality of floating-point numbers is fraught with peril. If the vector $u$ turns out to be very near to $-\zv$, then the remaining computations may be numerically a bit unstable, and you really don't want that. On the other hand, if you're just making a model of your back porch to decide where to put the picnic table, this nasty case is unlikely to affect you, and you can stop reading. 

To continue: to make this whole approach a little more numerically robust, I'll modify things from step 2 onwards: if the vector $u$ is too close to being $\pm \zv$, I'm going to rotate the cylinder to point along the $y$-axis first. I'll use $\yv$ to denote the vector $(0, 1, 0)$. Here's the revised version, math only. I leave it to you to translate to openSCAD. 


*

*Compute $v = Q - P$, the vector from the point $P$ to the point $Q$. If $v$ is zero, your input was bad (a length-zero cylinder makes no sense), and you'll have to handle this exceptional case. In fact, if $\| v \|$ is smaller than, say, $10^{-5}$, you should probably treat this as an error case. (I can't recall whether openSCAD uses IEEE floats or IEEE doubles or something else to represent things, but using $10^{-5}$ as a threshold should keep you completely safe regardless. You could probably get away with using $10^{-10}$.) 

*Compute $u = \frac{v}{ \| v\|}$. If the absolute value of the $z$-component of $u$ is no more than $0.8$, let $\sv = \zv$; if the absolute value of the $z$-component is more than $0.8$, then rotate the cylinder by $-90$ degrees about the $x$-axis, so that it now lies along the positive $y$ axis, and let $\sv = \yv$. In both cases, $\sv$ refers to "the starting orientation of the cylinder".

*Compute $w = \frac12 (\sv + u)$, and $t = \frac{w}{\|w\|}$,  At this point, $t$ is a unit vector in the plane of $v$ and $\sv$, and is exactly halfway between those two vectors. 

*Rotate about $t$ by 180 degrees to take the cylinder, which used to point along the $\sv$ direction, and make it point along the $v$ direction instead. 

*Scale everything by $d = \| v \| $, the distance from $P$ to $Q$. The cylinder now lies along the $PQ$ direction, with one end at the origin, and the other end $d$ units from the origin, i.e., it has the same length and orientation as the line segment $\overline{PQ}$. 

*Translate the cylinder by $P$ to move the origin (one end of the cylinder) to $P$, and the other end (which is at the tip of $v = Q - P$) to the point $Q$. 

There's another method that involves constructing a $4 \times 4$ matrix directly, and using openSCAD's multmatrix operation, but it's probably a little messier to write and to understand, so unless you're unhappy with the solution above, I won't bother writing it out. 
A: Start with your cylinder of radius $d$ and height $h.$
According to the OpenSCAD documentation, if you use the default value of center, that is, if it is false, you get a cylinder whose ends are circles with centers at $(0,0,0)$ and $(0,0,h)$.
Consider a vector between your two desired endpoints.
If you let $\vec p_1$ and $\vec p_2$ be the position vectors of the two endpoints,
then the vector you have defined is
$$\vec F = \vec p_2 - \vec p_1
= (x_2-x_1) \hat{\imath} +(y_2-y_1)\hat{\jmath} + (z_2-z_1)\hat{k}.$$
Let's find a rotation that turns the vector $h\hat k$ into $\vec F,$
that is, we want to rotate the vector that goes from the bottom to the top of the cylinder so that it points in the direction of $\vec F .$
One way to attack this is to work out $\vec F$ in spherical coordinates,
that is, find its angle of rotation $\phi$ away from the $z$ axis (its colatitude) and its angle of rotation $\theta$ around the $z$ axis (its longitude).
Mathematically, the colatitude is
$$ \phi = \arccos\left(\frac{z_2 - z_1}{\lVert \vec F \rVert} \right).$$
But if $\lvert z_2 - z_1\rvert$ is very close to $\lVert \vec F \rVert,$
the angle you get by this method is not very accurate.
There are various ways to work around this.
One way (based on the arctan formula on page 15 of this paper) is
$$
\phi = \mathrm{atan2} (\lVert \vec F - h\hat k\rVert, \lVert \vec F + h\hat k\rVert),
$$
where
\begin{align}
\lVert \vec F - h\hat k\rVert &= \sqrt{(x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1-h)^2},\\
\lVert \vec F + h\hat k\rVert &= \sqrt{(x_2-x_1)^2 + (y_2-y_1)^2 + (z_2-z_1+h)^2}.
\end{align}
and $\mathrm{atan2}$ is the usual two-argument arc tangent function.
(I omitted the step in Kahan's formula where each vector is divided by its magnitude;
since the magnitudes of the vectors are the same in this case, these factors cancel out anyway in the calculation of the arc tangent.)
For the longitude, you can take
$$
\theta = \mathrm{atan2}\left(\frac{y_2-y_1}{\sqrt{(x_2-x_1)^2 + (y_2-y_1)^2}},
\frac{x_2-x_1}{\sqrt{(x_2-x_1)^2 + (y_2-y_1)^2}}\right).
$$
Now to rotate the cylinder so it points in the right direction, you rotate by
$\phi$ around the $y$ axis and by $\theta$ around the $z$ axis.
Finally, translate by $\vec p_1$ to move the two ends of the cylinder to the desired points.
Note that since a cylinder is rotationally symmetric around its axis, you don't need three rotations around coordinate axes to orient it whichever way you want.
Two rotations is enough.

An alternative method is inspired by John Hughes' answer.
The idea is to rotate around a vector that is "midway" between $\vec F$ and the $z$ axis.
The vector $\vec F + h\hat k$ is one candidate, but it has an inaccuracy due to cancellation when the angle between the vectors is large.
So we make sure the angle between the vectors is not too large.
If $z_2 \geq z_1,$ then the angle between $\vec F$ and $h\hat k$ is a right angle or less. We can set
$$ \vec v =  h\hat k + \vec F, $$
rotate the cylinder by $\pi$ radians ($180$ degrees) around $\vec v,$
then translate by $\vec p_1$ so that the endpoints of the cylinder are where you want them.
If $z_2 < z_1$ then the angle between $\vec F$ and $h\hat k$ is greater than a right angle. We fix this by (in effect) reversing the direction of $\vec F$.
That is,  set
$$ \vec w =  h\hat k - \vec F, $$
and rotate the cylinder by $\pi$ radians  around $\vec w.$
This gives you a cylinder with one end at $(0,0,0),$ and the direction from there to the other end of the cylinder is the direction of $-\vec F,$ that is, the direction from $p_2$ to $p_1.$
The last step then is to translate by $\vec p_2$ (not $\vec p_1$) so that the end that is at $(0,0,0)$ moves to $p_2$ and the other end moves to $p_1.$
