Is this system of geometric points solvable? Apologies for the vague title, but I am not sure how otherwise to explain my question with a title. I have a set of 3 known points (x, y), and a set of 3 unknown points and a few known relationships between the points.

The positions of the red points are known.
The context is that a Quadratic Bezier curve drawn from s' to e' with control point a' will pass through a at its midpoint period (I've solved the math for this part already).

*

*s' is at the intersection of the line between a' and s and of the circle with known radius r around s.

*e' is at the intersection of the line between a' and e and of the circle with known radius r around e.

*a' can be solved from the three unknowns: $a_x = (2a'_x + s'_x + e'_x) / 4$. With the same for the y coordinate.

Putting that all together, I need a' to solve for s' and e', but I need those two terms to solve for a'. I've tried to put together a system of equations.

*

*$s'_x = s_x + r * cos(arctan(\frac{a'_y - s_y}{a'_x - s_x}))$

*$s'_y = s_y + r * sin(arctan(\frac{a'_y - s_y}{a'_x - s_x}))$

*$e'_x = e_x + r * cos(arctan(\frac{a'_y - e_y}{a'_x - e_x}))$

*$e'_y = e_y + r * sin(arctan(\frac{a'_y - e_y}{a'_x - e_x}))$

*$a'_x = (2*a_x + s'_x + e'_x) / 4$

*$a'_y = (2*a_x + s'_x + e'_x) / 4$
So I have a system of 6 equations with 6 unknowns that I am not able to solve. The trig gets way too crazy for me, and I am sure that I am missing many simplification steps or tricks somewhere.
In some situations, I know a' instead of a, and in those cases, I can solve for a. So I believe there should be at least (and probably only) one solution to the problem above that I can't solve. Also, there is a situation with doing a sign flip somewhere after the arctan that I've figured out, but not sure how to represent in math notation.
Any help is appreciated. Thanks in advance.

Supplementary information:
For even more context, in case anyone ever runs into this. I am
designing a graph (nodes and edges graph) visualization interface,
where I would like to allow the user to bend the edges manually when
the layout algorithm has not produced a desirable result (as is
common). The offsets s' and e' are  to account for the node radii
such that the terminal point of the edge is always directed to the
center of the node (which, especially with arrow-heads, looks the
best). a is fixed to the mouse position, but the svg engine (and underlying Bezier function of course) needs
to be specified the a' parameter. My current implementation was to
start off with a straight line where $a = a'$, and at each iteration
of edge-bending, to approximate a' by assuming the s' and e' of
the previous iteration. Basically, $s'_i = s'_{i-1}$ for the purposes of calculating $a'_i = (2 * a + s'_{i-1} + e'_{i-1}) / 4$. This performs pretty well in most cases, but
there are certain regions of a where the approximation becomes
wildly off and causes bouncing-back and forth effect until a moves
out of the region.
The region of terrible-approximation is something like this:
, and I would be OK with falling back to a
straight-line in this admittedly uninteresting region, but I also
can't seem to detect the boundaries of the region either.

 A: Treat all points as complex numbers. There is then really only one variable $a'$, for the curve endpoints may be derived from it as follows.
$s'$ lies on the line between $s$ and $a'$. The distances $ss'$ and $sa'$ are known (in terms of $a'$), so writing $u=\frac r{|s-a'|}$ as the ratio of those two distances, we have $s'=(1-u)s+ua'$. Similarly, if $v=\frac r{|e-a'|}$ then $e'=(1-v)e+va'$.
Then the single complex equation to solve is $s'+2a'+e'=4a$. Since the number of variables and equations is the same whether you split this into real and imaginary parts or not, the system is solvable… but in general you have to solve a sextic for each of real and imaginary parts, so you are better off solving numerically. An mpmath program doing this is below.

#!/usr/bin/env python3
from mpmath import *
mp.dps = 50

def solve(s0, a0, e0, r):
    def f(p): # p = (x, y)
        u = r / abs(p - s0)
        s1 = (1-u)*s0 + u*p
        v = r / abs(p - e0)
        e1 = (1-v)*e0 + v*p
        return s1 + 2*p + e1 - 4*a0
    # get a reasonable starting point
    u = r / abs(a0 - s0)
    s1 = (1-u)*s0 + u*a0
    v = r / abs(a0 - e0)
    e1 = (1-v)*e0 + v*a0
    p0 = (4*a0 - s1 - e1) / 2
    return findroot(f, p0, method="muller")

print(solve(-4+1j, 1+5j, 3-2j, 1))

