Given two points $p_1,p_2$ on a circular billard table. I want to know all billard trajectories from $p_1$ to $p_2$ hitting the boundary precisely once.
Model of the circular billard:
Denote the orgin of $\mathbb{R}^2$ by $M:=(0,0)$ and assume the billard table has radius one. Thus the boundary of the table is given by $\lbrace (x,y) \in \mathbb{R}^2 : x^2 +y^2 =1 \rbrace$. One can assume that the point $p_1$ is given by $p_1=(0,a)$ for some $a \in (-1,0)$. Notice that the problem is trivial if $p_1$ is either the midpoint $M$ or is lying on the boundary (i.e. $a=0$ or $a=-1$). The situation above can be achieved by rotation.
The point $p_2$ can after rotation and reflection assumed to be given by $p_2=(b,c)$ for some $(b,c) \in (0,1) \times (-1,1)$ satisfying $b^2 + c^2 \leq 1$.
Now the problem is to find for given $a,b,c$ a point $(x,y)$ such that the billard tratjectory from $p_1=(0,a)$ to $(x,y)$ is going through $p_2=(b,c)$. (Billard trajectory means that the angle of incidence is equal to the angle of reflection).
First approach:
Denote by $B$ the point at the boundary of the billard table such that the trajectory from $p_1$ to $p_2$ goes through $B$. For given $p_1=(0,a), p_2=(b,c)$ the distance from $p_1$ to $B$ and the distance from $B$ to $p_2$ should be minimal. Therefore to find $B$ define the function
$$f: S_1((0,0)) \rightarrow \mathbb{R}; x \mapsto \Vert p_1 - x \Vert + \Vert p_2 - x \Vert$$
Thus the problem above reduces to finding the minimum of $f$. This problem can be attacked using Lagrange multipliers.
Thus one can assume $f$ is defined an all of $\mathbb{R}^2$. Its gradient is then given by
$$ \nabla f(x,y) = \begin{pmatrix} \frac{x}{\sqrt{x^2 + (y-a)^2}} + \frac{x-b}{\sqrt{(x-b)^ 2+ (y-c)^2}}\\ \frac{y-a}{\sqrt{x^2 + (y-a)^2}} + \frac{y-c}{\sqrt{(x-b)^2 + (y-c)^2}} \end{pmatrix}. $$
At this point I get stuck, since calculations seem to get to difficult.
I prefer to choose coordinate axes such that points $p_1$ and $p_2$ have the same $y$: $$ p_1=(a,d),\quad p_2=(b,d). $$ Let $M=(0,0)$ and $V=(x,y)$ be the bouncing point ($x^2+y^2=1$), so that $VM$ is the bisector of angle $\angle p_1Vp_2$. If $N$ is the intersection between $VM$ and $p_1p_2$, then: $N=(dx/y,d)$. On the other hand, from the bisector theorem we have $Vp_1:Vp_2=Np_1:Np_2$, that is: $$ [(x-a)^2+(y-d)^2](by-dx)^2=[(x-b)^2+(y-d)^2](ay-dx)^2. $$ This equation has the spurious solution $y=d$: factoring that out and performing some simplification leads to: $$ 2d(a+b)y^2-2(d^2-ab)xy+2dx-(a+b)y-(a+b)d=0. $$ Together with $x^2+y^2=1$, this equation forms a fourth order system, which can have up to four real solutions. On the other hand, one can easily find some examples where four solutions actually exist, for instance: $a = -0.7$, $b = 0.6$, $d = 0.1$.
An explicit solution in terms of $a$, $b$, $d$ is too large to be of any use. You'd better substitute into the equation the values of the parameters before solving it.
EDIT.
The above equation represents a hyperbola and the bouncing points are then the intersections between this hyperbola and the unit circle. They can thus be found with GeoGebra: https://www.geogebra.org/m/WW3KkeJB.