Find the point equidistant from two points and a line

Assuming such point $Q$ exists it must lie on the Bisector Line b of $P_1$ and $P_2$ i.e. the line through the midpoint of $P_1$ and $P_2$ and orthogonal to the line $\vec {P_1P_2}$.

Thus you can write down a parametric expression for $Q=Q(s)\in b$ and set the following equation for distances:

$$d^2(Ql)=d^2(QP_1)$$

Firstly note that, without loss of generality, we can assume that the origin coincides with the intersection point between $l$ and $b$. In fact if $b\equiv l$ solution is trivial otherwise if not we can find the intersection point and simply shift the axes (we could also rotate the axes in such way that $l$ or $b$ coincide with an axis but it should be more complex whereas shifting is trivial).

Thus, let's assume:

$P_1=(x_1,y_1), P_2=(x_2,y_2)$,

l: $ ax+by=0$,

b: $cx+dy=0$

NOTE

find c and d is trivial

Finding perpendicular bisector of the line segement joining $ (-1,4)\;\text{and}\;(3,-2)$

Parametric equation of $Q \in b$ is given by:

$Q(t)=Q(d \cdot t,-c \cdot t)$ with $t\in \mathbb{R}$

The distances are given by:

$$\text{d($Ql$)} = \frac{\left | Ax_{0} + By_{0} + C\right |}{\sqrt{A^2 + B^2} }= \frac {\left| ad\cdot t - bc\cdot t \right|}{\sqrt{a^2 + b^2}} $$

$$\text{d($QP_1$)} = \sqrt{(d\cdot t-x_1)^2 + (-c\cdot t-y_1)^2}$$

and thus

$$\frac {\left| adt - bct \right|}{\sqrt{a^2 + b^2}}=\sqrt{(d\cdot t-x_1)^2 + (-c\cdot t-y_1)^2}$$

$$\frac {\left( ad\cdot t - bc\cdot t \right)^2}{{a^2 + b^2}}=(d\cdot t-x_1)^2 + (-c\cdot t-y_1)^2$$

from which "t" and thus "Q" can be easily found.

Without loss of generality, $l$ is the $x$ axis, $p_1$ is at $(-a,y_1)$ and $p_2$ at $(a,y_2)$.

One of the parabolas is $$y^2=(x+a)^2+(y-y_1)^2$$ and the other

$$y^2=(x-a)^2+(y-y_2)^2.$$

After elimination of $y$,

$$y_2(x+a)^2-y_1(x-a)^2+y_1y_2(y_1-y_2)=0$$ can give you two solutions.

An approach that could be useful not only to provide the solutions, but also to discuss and understand their existence, is as follows. Given two points $P_1$ and $P_2$ and a line, we can set, without loss of generality, an $xy$ plane where the $x$-axis coincides with the line and the non-negative portion of the $y$-axis contains $P_1$. The coordinates of $P_1$ and $P_2$ can then be written as $(0,y_1)$ and $(x_2,y_2)$, respectively (with $y_1 \geq 0\,$). The segment $P_1P_2$ has slope $(y_2-y_1)/x_2$ and its middle point has coordinates $(x_2/2,(y_1+y_2)/2)\,\,\,\,$. So the equation of its geometric axis is

$$y=-\frac{x_2}{y_2-y_1} x + \frac{y_1+y_2}{2}+\frac{x_2^2}{2(y_2-y_1)} $$

Now we have to identify, on this line, a point whose distance from $P_1$ (or $P_2$) is equal to the distance from the $x$-axis. Therefore, if we call $(X,Y)$ the coordinates of this point, we must have

$$X^2+(Y-y_1)^2=Y^2$$

Since $(X,Y)$ must also satisfy the equation of the geometric axis, solving the system and simplifying we get that the searched coordinates are

$$X=\frac{y_1 x_2}{y_1 - y_2} \pm \sqrt{y_1 y_2\left(1 + \frac{x_2^2}{(y_1 - y_2)^2} \right) }$$ $$ Y= \frac{y_1+y_2}{2} + \frac{x_2 (2X-x_2)}{2(y_1 - y_2)}$$

Note that these formulas are valid for $y_1 \neq y_2\,$. If $y_1=y_2\,$ (i.e., if the geometric axis is perpendicular to the $x$-axis), we trivially have

$$X=\frac{x_2}{2}$$ $$Y=\frac{y_1}{2}+\frac{X^2}{2y_1}$$

Also note that, for the solutions to be real, $y_1y_2$ must be $\geq 0\,\,$. Because $y_1$ has been assumed $ \geq 0$ in the initial assumptions ($P_1$ is on the non-negative portion of the $y$-axis) this implies that:

• in the case $y_1 \neq y_2 \,$, if $y_1$ and $y_2$ are both $>0$ (i.e., $P_1$ is on the positive $y$-axis and $P_2$ is in the first or second quadrant) we have two different real solutions;

• in the case $y_1 \neq y_2\,$, if $y_1=0$ or $y_2=0$ (i.e., one among $P_1$ and $P_2$ is on the $x$-axis) we have a single real solution;

• in the case $y_1 \neq y_2\,$, if $y_1>0$ and $y_2<0$ (i.e., $P_1$ is on the positive $y$-axis and $P_2$ is in the third or fourth quadrant) we have no real solutions;

• lastly, in the case $y_1 = y_2\,$, we have a single real solution.

To provide an example, let us set $P_1=(0,2) \,\,$ and $P_2=(2,4)\,\,$. In this case, we have $y_1=2\,$, $x_2=2\,\,$,$y_2=4\,\,$, and we are in the case $y_1 \neq y_2\,\,\,$. Applying the formulas above, we get

$$X=\frac{4}{-2} \pm \sqrt{2 \cdot 4 \left(1 + \frac{2^2}{(-2)^2} \right) }=-2 \pm 4$$

Putting these two possible values of $X$ in the equation giving $Y$, for the case $X=2\,\,$ yields

$$ Y= \frac{2+4}{2} + \frac{2 (4-2)}{2(-2)}= 3-1=2$$

and for the case $X=-6 \,\,$ yields

$$ Y= \frac{2+4}{2} + \frac{2 (2 \cdot (-6)-2)}{2(-2)}= 3+7=10$$

finally giving the two solutions $(2,2) \,$ and $(-6,10)\,$. Accordingly, the distances of $(2,2) \,$from $P_1$, $P_2$, and the $x$-axis are all equal to $2$, and those of $(-6,10)\,$ are all equal to $10$.

In the same manner, if for example we set $P_1=(0,4)\,\,$ and $P_2=(12,4)\,\,$, then $y_1=4\,$, $x_2=12\,\,\,$,$y_2=4\,\,$ and we are in the case $y_1 = y_2\,\,$, where the geometric axis is perpendicular to the $x$-axis. So we obtain

$$X=\frac{12}{2}=6$$ $$Y=\frac{4}{2}+\frac{6^2}{2 \cdot4}=\frac{13}{2}$$

giving the single solution $(6,13/2) \,\,$, whose distances from $P_1$, $P_2$, and the $x$-axis are all equal to $13/2 \,$.

I think you are almost there. The parabola has a property you already know.There are $two$ solutions/points for circle centers, but not one, get detected by a direct procedure as follows:

Intersections of a properly/conveniently placed parabola ( wlog $y=-f$ is chosen directrix) and perpendicular bisector of $P_1P_2 $ should be found.

$$ P_1:(a,b)\, ; P_2: (0,f) ; $$ where $f$ is parabola focal length.

Bisector equation is

$$(x-a)^2+(y-b)^2= x^2+(y-f)^2$$

Simplifying

$$ y(f-b)-ax+Q=0 ,\,where\, Q=(a^2+b^2-f^2)/2 \tag 1 $$

Parabola equation

$$ y = \frac{x^2}{4f} \tag2 $$

Plug 2) into 1) and solve quadratic in $x$, getting two solutions, meaning two points satisfy the given condition.

The key point to realize is... there is a single parabola, not two, containing centers of these two circles on the parabola, but not one circle.

$$ x_1 = 2/ (1 - b/f)*(a + \sqrt{a^2 - Q (1 - b/f)}; \, y_1 = x_1^2/(4 f) ; $$

$$ x_2 = 2/ (1 - b/f)*(a - \sqrt{a^2 - Q (1 - b/f)}; \, y_2 = x_2^2/(4 f); $$

For the numerical example given in graph, I have taken the values $ (a,b,f)=((1,3,1)$ as numerical input.

and it results in coordinates of intersection points of circumcircle centers. By evaluation/calculation approximately they are :

$$C_1=(-4.16228,4.33114);\, C_2=(2.16228,1.16886) ; \,$$ which depicts two circumcircles as you desired.

EDIT1:

Another possibility for $l$ seems to exist as given schematically; however it is not an independent situation for the line given but reflection about line of centers.

EDIT2:

If we take $(P_1,P_2)$ as $(-c,0),(c,0)$ and the given line as $ y= m x + c, $ the equations to the (tilted double parabola ) would be simplest.

I personally prefer to attack such problems first using basic vector algebra, mostly because they tend to be easier to implement for numerical computations.

In other words, this answer assumes the underlying problem to be solved, is to find and implement a stable numerical algorithm for solving specific problems of this type; and not something like a geometric exercise or educational investigation.

First, we want to simplify the situation. The obvious operation is to transform the system so that the line $l$ is on the $x$ axis.

Let's define the line $l$ in vector form: $$\vec{l}(t) = \vec{l}_0 + t \hat{l}_n \tag{1}\label{NA1}$$ where $\vec{l}_0 = (x_0 ,\, y_0)$ is any point the line passes through, and $\hat{l}_n = (x_n ,\, y_n)$ is an unit vector ($x_n^2 + y_n^2 = 1$).

If we define rotation matrix $\mathbf{R}$, $$\mathbf{R} = \left [ \begin{array}{cc} x_n & y_n \\ y_n & -x_n \end{array} \right ] \tag{2}\label{NA2}$$ then the transformation we want is $$\vec{p}' = \mathbf{R}\left(\vec{p} - \vec{l}_0\right)\tag{3a}\label{NA3a}$$ i.e. point $\vec{p} = (x, y)$ becomes $$\left\lbrace\begin{array}{l} x' = x_n ( x - x_0 ) + y_n ( y - y_0 ) \\ y' = y_n ( x - x_0 ) - x_n ( y - y_0 ) \end{array}\right.\tag{3b}\label{NA3b}$$

The inverse operation is $$\vec{p} = \vec{l}_0 + \mathbf{R}^T \vec{p}'$$ and because $\mathbf{R}^T = \mathbf{R}$ in this particular case, we can write it as $$\left\lbrace\begin{array}{l} x = x_0 + x_n x' + y_n y' \\ y = y_0 + y_n x' - x_n y' \end{array}\right.\tag{3c}\label{NA3c}$$

Next, we can consider this much simpler situation.

Let's say the two points (as transformed by $\eqref{NA3b}$) are $\vec{p}_1 = ( x_1 ,\, y_1 )$ and $\vec{p}_2 = ( x_2 ,\, y_2 )$.

The requirement for the point $\vec{p} = ( x ,\, y )$ we seek is that it is as far from $\vec{p}_1$ and $\vec{p}_2$ as it is from the line $l$, i.e. from $y = 0$. $$\begin{cases} \lvert y \rvert = \sqrt{ (x - x_1)^2 + (y - y_1)^2 } \\ \lvert y \rvert = \sqrt{ (x - x_2)^2 + (y - y_2)^2 } \end{cases}$$ Because both sides are nonnegative, we can square them, to get the pair of equations in a more tractable form, $$\begin{cases} y^2 = (x - x_1)^2 + (y - y_1)^2 \\ y^2 = (x - x_2)^2 + (y - y_2)^2 \end{cases} \iff \begin{cases} (x - x_1)^2 + y_1^2 - 2 y y_1 = 0 \\ (x - x_2)^2 + y_2^2 - 2 y y_2 = 0 \end{cases} \tag{4}\label{NA4}$$ After finding the solution or solutions $(x, y)$, use $\eqref{NA3c}$ to transform them back to the original coordinate system.

Let's consider the numerical solutions to $\eqref{NA4}$, and the cases where there is no solution at all. (I shall quietly assume you are aware of the pitfalls when using floating-point numbers; especially allowing for rounding error, "machine epsilon", when comparing values for equality.)

1. If $y_1 y_2 \lt 0$, there can be no result, because the given points are on different sides of the line: the distance from the tentative result point to at least one of the given points is always greater than the distance from the tentative result point to the line.

2. If $y_1 \ne y_2$, but $y_1 = 0$ or $y_2 = 0$ (or $y_1 y_2 = 0$), there can be no result, because only one of the given points is on the line $l$, and the common distance cannot be zero and nonzero at the same time.

3. If $y_1 = y_2 = 0$ but $x_1 \ne x_2$, there is no solution, because the given points are separate but both on the line $l$, and the common distance cannot be both zero and nonzero at the same time.

4. If $y_1 = y_2 = 0$ and $x_1 = x_2$, the given points are at the same point on the line $l$, so the common distance is zero. The result is the position of the given points: $$\begin{cases}x = x_1 = x_2 \\ y = y_1 = y_2 = 0 \end{cases}$$

5. If $y_1 = y_2 \ne 0$ and $x_1 = x_2$, the given points are the same, and the result is at midway between the given points and the line: $$\begin{cases}x = x_1 = x_2 \\ y = \frac{y_1}{2} = \frac{y_2}{2} \end{cases}$$

6. If $y_1 = y_2 \ne 0$ and $x_1 \ne x_2$, there is a single solution: $$\begin{cases}x = \frac{x_1 + x_2}{2} \\ y = \frac{y_1}{2} + \frac{(x_2 - x_1)^2}{8 y_1} \end{cases}$$

7. If $y_1 \ne y_2$, $y_1 y_2 \gt 0$, and $x_1 = x_2$, there given points are on the line perpendicular to line $l$, on the same side of line $l$, and there is a pair of solutions: $$\begin{cases} x = x_1 \pm \sqrt{y_1 y_2} \\ y = \frac{y_1 + y_2}{2} \end{cases}$$

8. Otherwise, $y_1 \ne y_2$, $y_1 y_2 \gt 0$, and $x_1 \ne x_2$, and there are two solutions: $$\begin{cases} x = \frac{x_1 y_2 - x_2 y_1 + \sqrt{y_1 y_2 \left ( (x_2 - x_1)^2 + (y_2 - y_1)^2 \right )}}{y_2 - y_1} \\ y = \frac{(y_1 + y_2) \left ( (x_2 - x_1)^2 + (y_2 - y_1)^2 \right ) - (x_2 - x_1)\sqrt{y_1 y_2 \left ( (x_2 - x_1)^2 + (y_2 - y_1)^2 \right )}}{(y_2 - y_1)^2 } \end{cases}$$ and $$\begin{cases} x = \frac{x_1 y_2 - x_2 y_1 - \sqrt{y_1 y_2 \left ( (x_2 - x_1)^2 + (y_2 - y_1)^2 \right )}}{y_2 - y_1} \\ y = \frac{(y_1 + y_2) \left ( (x_2 - x_1)^2 + (y_2 - y_1)^2 \right ) + (x_2 - x_1)\sqrt{y_1 y_2 \left ( (x_2 - x_1)^2 + (y_2 - y_1)^2 \right )}}{(y_2 - y_1)^2 } \end{cases}$$ If you need to choose only one of the solutions, I recommend choosing the one with minimum $y$ in magnitude, $\lvert y \rvert$, as that solution corresponds to the smaller common distance.

If there is a solution, you can use $\eqref{NA3c}$ to convert back to the original coordinate system.

Note that the above cases are not in recommended implementation order; if you draw the cases into a diagram, you'll see you have a simple decision tree. As conditionals tend to be slow on current hardware, implementing a decision tree (to find the solution in fewest average and maximum conditional tests) is recommended.

Also note that I obtained the above solutions using free and open source Sagemath, but there may be transcription errors. If you find an error, or a suspicious point, please let me know in a comment so I can verify and fix.