Orthogonal line segment between line segments with length I am given two line segments $\overline{AB}$ and $\overline{CD}$. I am looking for the line segment $\overline{EF}$ between $\overline{AB}$ and $\overline{CD}$ which has a desired length $l$ while also being orthogonal to $\overline{AB}$, whereas $E$ should lie on $\overline{AB}$ and $F$ on $\overline{CD}$. All of this is happening in $\mathbb{R^3}$. When multiple solutions exist, we want to know all of them. We also want to know when none exist.
To sum up:
$$(\overrightarrow{b}-\overrightarrow{a})\bullet(\overrightarrow{f}-\overrightarrow{e})=0$$
$$|\overrightarrow{f}-\overrightarrow{e}| = l \quad,l\ge0$$
$$\overrightarrow{e} = \overrightarrow{a}+m(\overrightarrow{b}-\overrightarrow{a}) \quad,0\le m\le1$$
$$\overrightarrow{f} = \overrightarrow{c}+n(\overrightarrow{d}-\overrightarrow{c}) \quad,0\le n\le1$$
We can insert the third and fourth requirement into the first to get:
$$(\overrightarrow{b}-\overrightarrow{a})\bullet(\overrightarrow{c}+n(\overrightarrow{d}-\overrightarrow{c})-\overrightarrow{a}-m(\overrightarrow{b}-\overrightarrow{a}))=0$$
From here on, however, I struggle to incorporate the second requirement. How do I get the length constraint in there as well? Here is a different approach which also leaves me stranded:
The wanted line segment $\overline{EF}$ lies somewhere on the plane with normal vector $\frac{\overrightarrow{b}-\overrightarrow{a}}{|\overrightarrow{b}-\overrightarrow{a}|}$ and origin somewhere on $\overline{AB}$. This yields the plane equation
$$\frac{\overrightarrow{b}-\overrightarrow{a}}{|\overrightarrow{b}-\overrightarrow{a}|}\bullet\overrightarrow{p}+d=0$$
where $\overrightarrow{p}$ is a point on the plane and $d$ is the (signed) distance of the plane from the origin.
For $\overrightarrow{e}$, this would be:
$$\frac{\overrightarrow{b}-\overrightarrow{a}}{|\overrightarrow{b}-\overrightarrow{a}|}\bullet(\overrightarrow{a}+m(\overrightarrow{b}-\overrightarrow{a}))+d=0$$
Since $\overrightarrow{e}$ lies on $\overline{AB}$, which is parallel to our plane normal, we can express it using $d$ instead of $m$, effectively getting rid of one unknown variable (mind the sign):
$$\frac{\overrightarrow{b}-\overrightarrow{a}}{|\overrightarrow{b}-\overrightarrow{a}|}\bullet(-d\frac{\overrightarrow{b}-\overrightarrow{a}}{|\overrightarrow{b}-\overrightarrow{a}|})+d=0$$
For $\overrightarrow{f}$ the plane equation is:
$$\frac{\overrightarrow{b}-\overrightarrow{a}}{|\overrightarrow{b}-\overrightarrow{a}|}\bullet(\overrightarrow{c}+n(\overrightarrow{d}-\overrightarrow{c}))+d=0$$
Here, we want to solve for $n$ to get rid of that unknown as well. This will allow us to define $\overrightarrow{e}$ and $\overrightarrow{f}$ using a single variable in such a way that $\overrightarrow{f}-\overrightarrow{e}$ will always express a vector which is orthogonal to the plane normal. To do this, we need to explode the equation:
$$\frac{(b_x-a_x)(c_x+n(d_x-c_x))+(b_y-a_y)(c_y+n(d_y-c_y))+(b_z-a_z)(c_z+n(d_z-c_z))}{\sqrt{(b_x-a_x)^2+(b_y-a_y)^2+(b_z-a_z)^2}}+d=0$$
Here, I struggle to solve for $n$.
If I didn't, though, I could then solve the following for $d$
$$l=|\overrightarrow{f}-(-d\frac{\overrightarrow{b}-\overrightarrow{a}}{|\overrightarrow{b}-\overrightarrow{a}|})|$$
where $\overrightarrow{f}$ would be replaced by an equation which expresses $\overrightarrow{f}$ through $d$ instead of $n$. I could then solve $\overrightarrow{e}$:
$$\overrightarrow{e}=-d\frac{\overrightarrow{b}-\overrightarrow{a}}{|\overrightarrow{b}-\overrightarrow{a}|}$$
and – if I had figured it out – $\overrightarrow{f}$.
 A: What you’re essentially trying to find are the intersections of $\overline{CD}$ with a piece of a cylinder of radius $l$ and axis $AB$. There can be zero, one, two or an infinite number of intersection points. I would suggest breaking this down into two stages. First, find all of the points on the lines $AB$ and $CD$ that meet the criteria—that is, the intersection of the line $CD$ with the entire cylinder—and then cull those solutions that don’t lie on the corresponding line segments.  
There are several ways to approach the unconstrained problem. A straightforward way is to parameterize $CD$ and then use a point-line distance formula. So, let $P(\lambda)=(1-\lambda)C+\lambda D$. From the distance formula we then get the equation $$\|(B-A)\times(P(\lambda)-A)\|^2 = l^2\|B-A\|^2. \tag 1$$ This will either be a quadratic equation in $\lambda$, or, if $CD$ is parallel to $AB$, an equality of two constants. Solutions in the interval $[0,1]$ lie on the segment $\overline{CD}$. The corresponding point on $AB$ is the intersection of this line with its perpendicular plane through $P(\lambda)$, which has equation $$(B-A)\cdot (X - P(\lambda))=0.\tag 2$$ This intersection point can be computed directly in various ways, such as the Plücker matrix of $AB$, but since you need to do a range check, anyway, parameterize $\overline{AB}$ as $Q(\mu)=(1-\mu)A+\mu B$ and substitute into (2) to get $$(B-A)\cdot Q(\mu)=(B-A)\cdot P(\lambda). \tag 3$$ This is a linear equation in $\mu$ that you can solve and again check that the solution lies in $[0,1]$. Alternatively, compare the signs of the values that you get by plugging $A$ and $B$ into the left-hand side of (2): if either is zero or they have opposite signs, then the intersection lies on the segment $\overline{AB}$ and you can go ahead and compute it using any convenient method.  
You could instead have started with equation (2), which expresses the constraint that the two points both lie on a plane perpendicular to $AB$. This is a linear equation in $\lambda$ and $\mu$. The distance constraint can be expressed as the equation $$\|P(\lambda)-Q(\mu)\|^2=l^2, \tag 4$$ so the problem is reduced to the two-dimensional problem of finding the intersections of a conic and line. There are algorithms for doing this that are suitable for automation if you don’t have access to an equation solver. Once again, after you’ve found the solutions to this system of equations, the ones you want are those for which $\lambda,\mu\in[0,1]$. As a third approach, you could derive the equation of the cylinder and set up equations for its intersection with $CD$, but that’s eventually going to lead to equations similar to those above with a lot more work.  
A: As amd suggested, a good starting point is the equation of an axis-aligned cylinder:
$$r^2=x^2+y^2$$
I ended up solving the problem by transforming $\overline{AB}$ and $\overline{CD}$ so that $A$ is at $(0,0,0)$ and $\overrightarrow{b}-\overrightarrow{a}$ points into positive z. I could then insert the $x$ and $y$ components from the definitions of $\overrightarrow{f}$ from my question into above formula and solve for $n$, determine transformed $\overrightarrow{f}$ and then transform the result back. I also defined:
$$\overrightarrow{v}=\overrightarrow{d}-\overrightarrow{c}$$
We then start with this:
$$l^2=(c_x+nv_x)^2+(c_y+nv_y)^2$$
Here, $r^2$ became $l^2$ from my question. Solving for $n$ using solver of choice yields:
$$n=\frac{\pm\sqrt{-c_x^2v_y^2+2c_xc_yv_xv_y-c_y^2v_x^2+l^2(v_x^2+v_y^2)}-c_xv_x-c_yv_y}{v_x^2+v_y^2}$$
A computational difficulty occurs when $\overline{AB}$ and $\overline{CD}$ are parallel, as it causes the divisor to approach zero. In or near these cases a distance test with any point on $\overline{CD}$ can determine whether it is at the wanted distance. Then we can determine the transformed $\overrightarrow{f}$, and transform it back.
Determining $\overrightarrow{e}$ is then as simple as projecting $\overrightarrow{f}$ onto $\overline{AB}$.
Finally, one needs to test that $\overrightarrow{e}$ lies on $\overline{AB}$ to validate the result.
