OK, to follow up on my last post, here's my progress.
Let's assume $\vec T_i$ are circle centers, $\vec n_i$ vectors normal to each circle, and $r_i$ their respective radii.
The first step is to project one circle onto the plane of another (in my case, circle 1 onto the plane of the circle 3). For sake of brevity I will use following function $project(\vec v_1,\vec v_2)=\vec v_1-(\vec v_2 \cdot \vec v_1) \cdot \vec v_2$.
$\vec T'_1 = project(\vec T_1,\vec n_3)$ - projection point of the center of the circle on the plane
$\vec T'_1 = \vec T_1 + \left[ \vec n_3 \cdot \left(\vec T_3 - \vec T_1 \right) \right] \cdot \vec n_3$ - projection point of the center of the circle on the plane
$\vec y_1 = \frac {project(\vec n_1,\vec n_3)}{\vert project(\vec n_1,\vec n_3) \vert}$ - y axis of the new coordinate system
$\vec x_1 = \frac {\vec n_3 \times \vec y_1}{\vert \vec n_3 \times \vec y_1 \vert}$ - x axis of the new coordinate system
$\vec x_1 = \frac {\vec y_1 \times \vec n_3}{\vert \vec y_1 \times \vec n_3 \vert}$ - x axis of the new coordinate system
$f = \vec n_1 \cdot \vec n_3$ - ratio of the minor axis to the radius of the circle (major axis)
Obtaining coordinates of the third circle in new (2D) coordinate system:
$ x_3 = \left( \vec T_3-\vec T'_1 \right) \cdot \vec x_1$
$ y_3 = \left( \vec T_3-\vec T'_1 \right) \cdot \vec y_1$
Now we have two parallel coordinate systems $x,y$ and $x_k,y_k$ with shift of $(x_3,y_3)$.
In $x,y$ system we have ellipse $\frac{x^2}{r_1^2} + \frac{y^2}{f^2r_1^2}=1$, which gives:
$y=\pm f\sqrt{r_1^2-x^2}$ , $\dot y= \mp \frac{f \cdot x}{\sqrt{r_1^2-x^2}}$ and $l=y-x \cdot \dot y=\pm \frac{f \cdot r_1^2}{\sqrt{r_1^2-x^2}}$
To see what would be l in the other coordinate system we 'calculate' $l'_k$ as:
$l'_k= x_3 \cdot \dot y + l - y_3 = \ldots = \mp f \frac {x_3 \cdot x - r_1^2}{\sqrt{r_1^2-x^2}}-y_3$
Doing the same for the circle $x_k^2+y_k^2=r_3^2$ in $x_k,y_k$ coordinate system results in following:
$y_k=\pm\sqrt{r_3^2-x_k^2}$, $\dot y_k=\mp \frac{x_k}{\sqrt{r_3^2-x_k^2}}$ and $l_k=y_k-x_k \cdot \dot y_k=\pm \frac{r_3^2}{\sqrt{r_3^2-x_k^2}}$
And now for the fun part- equalling the slopes of the tangents ($\dot y = \dot y_k$) will rid us of one unknown $x_k$:
$\frac{f \cdot x}{\sqrt{r_1^2-x^2}}=\frac{x_k}{\sqrt{r_3^2-x_k^2}}$
$\ldots$
$x_k = \pm \frac {f \cdot r_3 \cdot x}{\sqrt {r_1^2-\left(1-f^2\right)\cdot x^2}}$ (note: plus sign is for outer, minus for inner tangents)
Equating the $l'_k=l_k$ will give the 'final' equation to solve for $x$:
$\mp f \frac {x_3 \cdot x - r_1^2}{\sqrt{r_1^2-x^2}}-y_3=\pm \frac{r_3^2}{\sqrt{r_3^2-x_k^2}}$
the troublesome bit is getting rid of $\sqrt{r_3^2-x_k^2}$ which is:
$\sqrt{r_3^2-x_k^2} = r_3\frac{\sqrt{r_1^2-x^2}}{\sqrt {r_1^2-\left(1-f^2\right)\cdot x^2}}$
and with that we get the equation:
$\pm\left(x_3 \cdot x - r_1^2\right)-y_3\sqrt{r_1^2-x^2}=\pm r_3\sqrt {r_1^2-\left(1-f^2\right)\cdot x^2}$
Once we square that equation, sort it out a bit, and then square it again and sort out once again, we're left with:
$A x^4+B x^3 + C x^2 +D x + E = 0$
$A=\left[ x_3^2- y_3^2 + \left(1-f^2\right) r_3^2 \right]^2 + 4f^2 x_3^2 y_3^2$
$B=-4f^2 r_1^2 x_3 \left[ x_3^2 + y_3^2 + \left(1-f^2\right) r_3^2 \right]$
$C= 2 r_1^2 \left( f^2 r_1^2 +y_3^2 - r_3^2 \right) \left[ x_3^2- y_3^2 + \left(1-f^2\right) r_3^2 \right] + 4f^2 r_1^2 \left[ r_1^2 \left( f^2 x_3^2 + y_3^2 \right) - x_3^2 y_3^2 \right]$
$D= 4f^2 r_1^4 x_3 \left( -f^2 r_1^2 +y_3^2 + r_3^2 \right)$
$E= r_1^4 \left[ \left( f^2 r_1^2 +y_3^2 - r_3^2 \right)^2 -4 f^2 r_1^2 y_3^2 \right]$
Solving this equation will give four values of x on the projected ellipse (y coordinates should be calculated, bearing in mind that, generally speaking, $x_3 \cdot y_3>0 \Rightarrow x \cdot y<0$, and vice versa, while $x_k$ and $y_k$ should also be calculated, $x \cdot x_k >0$ for outer, and $x \cdot x_k<0$ for innner tangents).
This method so far has been verified with CAD models. What is left to do is to 'unproject' points on the ellipse to the original plane (plane of circle 1) and with corresponding points on the circle 3 we have four lines, and each of those lines will lie in two of the eight planes. To solve the rest of the problem, I think the best course of action would be to pass a plane perpendicular to each line through third circle (circle 2) centerpoint, project that circle onto that plane and find tangents to that ellipse from the point in which line intersects the plane. Getting the equations of planes should be straightforward from there... If anyone's still interested, I'll post the following procedure when I finish and check it.
EDIT I have corrected first and third equation because, even though they worked in 2D, this is the correct form and works in 3D. Here's the rest of procedure (including solving the equation)...
$p_1 = 2C^3 - 9BCD + 27AD^2 + 27BE^2 - 72ACE$
$p_2 = p_1 + \sqrt {-4 \left( C^2 - 3BD + 12AE \right) ^3 + p_1^2}$
$p_3 = \frac {1}{3A} \left( 2 \frac { C^2 - 3BD + 12AE}{p_2} \sqrt[3] {\frac {p_1}{2}} +1 \right) \sqrt[3] {\frac {p_1}{2}}$
$p_4 = \sqrt {\frac {B^2}{4A^2} - \frac{2C}{3A} + p_3}$
$p_5 = \frac {B^2}{2A^2} - \frac{4C}{3A} - p_3$
$p_6 = \frac {-B^3 + 4ABC - 8A^2D}{4A^3p_4}$
$_1xT_2 = -\frac {B}{4A} - \frac {p_4}{2} \pm \frac {\sqrt {p_5 - p_6}}{2}$ $\quad , \quad$
$_3xT_4 = -\frac {B}{4A} + \frac {p_4}{2} \pm \frac {\sqrt {p_5 + p_6}}{2}$
For each od those $xT_i$ where $i=1..4$, we calculate:
$yT_i=\pm f\sqrt{r_1^2-xT_i^2}$ $\quad , \quad$ ${x_k}_i = \pm \frac {f \cdot r_3 \cdot xT_i}{\sqrt {r_1^2-\left(1-f^2\right)\cdot xT_i^2}}$ $\quad , \quad$ ${y_k}_i=\pm\sqrt{r_3^2-{x_k}_i^2}$
$\vec {{TP_3}_i} = \vec T_3 + {x_k}_i \cdot \vec x_1 +{y_k}_i \cdot \vec y_1$ - tangent point on the circle 3
$\vec {{TP´_1}_i} = \vec T´_1 + xT_i \cdot \vec x_1 +yT_i \cdot \vec y_1$ - tangent point on projected ellipse
$\vec {{TP_1}_i} = \vec {{TP´_1}_i} + \frac {\vec n_1 \left( \vec T_1 - \vec {{TP´_1}_i} \right) }{\vec n_1 \cdot \vec n_3} \vec n_3$ - tangent point on circle 1 (point on the ellipse 'unprojected' onto the plane of circle 1)
$\vec {{l_1}_i} = \frac {\vec {TP_3}_i - \vec {TP_1}_i}{\vert \vec {TP_3}_i - \vec {TP_1}_i \vert}$ - tangent line on circles 1 and 3
Now we set up a plane (normal to line $\vec {{l_1}_i}$, through point $T_2$) and coordinate system:
$\vec {{y_2}_i} = \frac {project(\vec n_2, \vec {{l_1}_i})}{\vert project(\vec n_2, \vec {{l_1}_i}) \vert}$ $\quad , \quad$ $\vec {{x_2}_i} = \frac {\vec {{y_2}_i} \times \vec n_2}{\vert \vec {{y_2}_i} \times \vec n_2 \vert}$
$\vec {{TP´_3}_i} = \vec {{TP_3}_i} + \left[ \vec {{l_1}_i} \cdot \left( \vec T_2 - \vec {{TP_3}_i} \right) \right] \cdot \vec {{l_1}_i}$ - projection of $\vec {{TP_3}_i}$ onto the plane
Coordinates of projected point in 2D coordnate system:
${x_{TP}}_i = \left( \vec {{TP´_3}_i} - \vec T_2 \right) \cdot \vec {x_2}_i $
${y_{TP}}_i = \left( \vec {{TP´_3}_i} - \vec T_2 \right) \cdot \vec {y_2}_i $
$f´_i = \vert \vec {{l_1}_i} \cdot \vec n_2 \vert$ - ratio of the minor axis to the radius of the circle (major axis)
Now we can calculate coordinates of tangent points (rays from ${TP´_3}_i$) on the ellipse:
$_1{{xT2}_i}_2 = \frac {f´_i^2 r_2^2 {x_{TP}}_i - r_2 {y_{TP}}_i \pm \sqrt {f´_i^2 \left( {x_{TP}}_i^2 - r_2^2\right) + {y_{TP}}_i^2}}{f´_i^2 {x_{TP}}_i^2 + {y_{TP}}_i^2}$
For each of those $xT2_i$ we get corresponding $yT2_i = \pm f´_i \sqrt{r_2^2 - xT2_i^2}$. And each of those points will be 'unprojected' onto the circle plane.
$\vec {TP´_2}_i = \vec T_2 + xT2_i \cdot \vec {x_2}_i + yT2_i \cdot \vec {y_2}_i$
$\vec {{TP_2}_i} = \vec {TP´_2}_i + \left[ \vec n_2 \cdot \left( \vec T_2 - \vec {TP´_2}_i \right) \right] \cdot \vec n_2$
$\vec {l_2}_i = \frac {\vec {{TP_3}_i} - \vec {{TP_2}_i}}{\vert \vec {{TP_3}_i} - \vec {{TP_2}_i} \vert}$ - second line in the tangent plane (sic. there are two of those lines for each i... resulting in eight planes)
$\vec {n_1}_i = \frac {\vec {l_1}_i \times \vec {l_2}_i}{\vert \vec {l_1}_i \times \vec {l_2}_i \vert}$ - normal vector to the plane
${D_1}_i = -\vec {n_1}_i \cdot \vec {{TP_3}_i}$ - plane constant
$\vec {n_1}_i \cdot \vec T + {D_1}_i = 0$ - equation of tangent planes