An identity that might prove useful in this problem is
$$
\cot\left(\frac{\theta}{2}\right)=\frac{\sin(\theta)}{1-\cos(\theta)}=\frac{1+\cos(\theta)}{\sin(\theta)}\tag{1}
$$
In $\mathbb{R}^3$, one usually uses the cross product to compute the sine of the angle between two vectors. However, one can use a two-dimensional analog of the cross product to do the same thing in $\mathbb{R}^2$.
$\hspace{5cm}$
In the diagram above, $(x,y)\perp(y,-x)$ and so the $\color{#FF0000}{\text{red angle}}$ is complementary to the $\color{#00A000}{\text{green angle}}$. Thus,
$$
\begin{align}
\sin(\color{#FF0000}{\text{red angle}})
&=\cos(\color{#00A000}{\text{green angle}})\\[6pt]
&=\frac{(u,v)\cdot(y,-x)}{|(u,v)||(y,-x)|}\\[6pt]
&=\frac{uy-vx}{|(u,v)||(x,y)|}\tag{2}
\end{align}
$$
$uy-vx$ is the normal component of $(u,v,0)\times(x,y,0)$; thus, it is a two dimensional analog of the cross product, and we will denote it as $(u,v)\times(x,y)=uy-vx$.
Let $u_a=\frac{Q-P}{|Q-P|}$ and $u_b=\frac{R-P}{|R-P|}$, then since
$$
|A-P|=|B-P|=r\cot\left(\frac{\theta}{2}\right)
$$
we get
$$
\begin{align}
A
&=P+ru_a\cot\left(\frac{\theta}{2}\right)\\
&=P+ru_a\frac{1+u_a\cdot u_b}{|u_a\times u_b|}
\end{align}
$$
and
$$
\begin{align}
B
&=P+ru_b\cot\left(\frac{\theta}{2}\right)\\
&=P+ru_b\frac{1+u_a\cdot u_b}{|u_a\times u_b|}
\end{align}
$$
where we take the absolute value of $u_a\times u_b$ so that the circle is in the direction of $u_a$ and $u_b$.
Here's a completely different approach that might be easier to encode.
Your two given points ($(x_1, y_1)$ and $(x_2, y_2)$) and the
centers of the two desired circles are at the four vertices of a
rhombus with side length $r$.
You can use the Pythagorean Theorem to find the length of the diagonal
of the rhombus from $(x_1, y_1)$ to $(x_2, y_2)$.
Better still, divide by two so you now have half the length of the diagonal.
The two diagonals of a rhombus are perpendicular,
so the point $(x_1, y_1)$, the center of the rhombus, and one of the
circles' centers make a right triangle with hypotenuse $r$
and one leg equal to half the known diagonal.
Use the Pythagorean Theorem to find half the length of the other diagonal.
Now you just need to construct line segments of that length with one
end at the center of the rhombus, perpendicular to the known diagonal,
and the other end of each segment will be the center of one of the
desired circles.
That's the entire rationale of the procedure.
For the detailed calculations, I'll assign names to various lengths as
we go along in order to keep the equations from getting too ugly.
(You'd probably want to do this anyway when you do this in software.)
Let $x_a = \frac12(x_2 - x_1)$ and $y_a = \frac12(y_2 - y_1)$;
then the center of the rhombus is at $(x_0,y_0) = (x_1+x_a, y_1+y_a)$
and half the length of the known diagonal is
$a = \sqrt{x_a^2 + y_a^2}.$
The length of the other diagonal is
$b = \sqrt{r^2 - a^2}.$
Now suppose (for example) that $(x_0,y_0)$ is above and to the right
of $(x_1,y_1)$ (assuming, at least for now,
the usual mathematical Cartesian $x,y$ coordinates).
To get from $(x_1,y_1)$ to $(x_0,y_0)$ along a straight line,
we go to the right $x_a$ units and up $y_a$ units for every
$a$ units of movement.
So to get from $(x_0,y_0)$ to the center of one of the circles,
$(x_3,y_3)$, we can just turn those rules $90$ degrees:
we can go down $x_a$ units and to the right $y_a$ units for every
$a$ units of movement;
that is, $\frac{x_a}{a}$ down and $\frac{y_a}{a}$ right for each unit
traveled along the diagonal.
But the distance from $(x_0,y_0)$ to $(x_3,y_3)$ is $b$,
so we will end up going down $b\frac{x_a}{a}$ units
and to the right $b\frac{y_a}{a}$ units.
In other words,
\begin{align}
x_3 &= x_0 + \frac{b y_a}{a}, \\
y_3 &= y_0 - \frac{b x_a}{a}. \\
\end{align}
To get to the center of the other circle, $(x_4,y_4)$,
we just go the same amount in the opposite direction from $(x_0,y_0)$:
\begin{align}
x_4 &= x_0 - \frac{b y_a}{a}, \\
y_4 &= y_0 + \frac{b x_a}{a}. \\
\end{align}
If $(x_0,y_0)$ is not above and to the right of $(x_1,y_1)$,
or if you're working in a graphics coordinate system with $y$ downward,
or both, the signs of either $x_a$ or $y_a$ (or both)
might be negative instead of positive; but all the formulas above
continue to work exactly the same way as before.
Best Answer
Let's say $\epsilon$ is the line defined by A and B and $K_1,K_2$ and $K_3$ the centers of the given circles.
Then we have the following 3 equations that will help us define $x_3,y_3,r_3$
$1$. $\sqrt{(x_3-x_2)^2+(y_3-y_2)^2}=r_3+r_2$ or $\sqrt{(x_3-x_2)^2+(y_3-y_2)^2}=|r_3-r_2|$
That is, the distance of the centers $K_2,K_3$ equals the sum of the radius $r_3,r_2$ if the circles are tangent outwardly or the $|r_3-r_2|$ if they're tangent inwardly.
$2$. $\sqrt{(x_3-x_1)^2+(y_3-y_1)^2}=r_3+r_1$ or $\sqrt{(x_3-x_1)^2+(y_3-y_1)^2}=|r_3-r_1|$
$3$. The distance of the $K_3$ from $\epsilon$ equals $r_3$ if you write $\epsilon :Ax+By+C=0$ then $r_3=\frac{|Ax_3+By_3+c|}{\sqrt{A^2+B^2}}$