Find point of reflection on circle I'm playing around with some ray-tracing type applications and I've run into following problem:

I'm given two points $A$ and $B$ and a circle with the center at $(0,0)$ and radius $r$. $A$ and $B$ have rays that meet and touch the circle at $P$. How can I find P such that the angles of incidence for the two rays are the same (ie $\alpha$ = $\beta$)? $A$ and $B$ are guaranteed to be 'nice', ie they aren't on opposite sides of the circle and the angle of incidence is greater than 0.
Currently I am solving this numerically, but I feel sure that there must be a closed form solution to this, I'm just having some trouble finding it.
 A: Here is a potential approach.
By rotation, you can assume that $A = (0, h)$ for some $h>r$, and that $B = (x, y)$ for some $x>0$ (and $x^2 + y^2 > r$).
Now, let $P = (a, b)$ be a point on the circle (with $a>0$).
Then the angle $\alpha$ should be formed by the vectors $\langle -b, a \rangle$ and $\langle -a, h-b \rangle$, whose dot product one can readily compute. And the angle $\beta$ should be formed by $\langle b, -a \rangle$ and $\langle x-a, y-b \rangle$, again whose dot product can be computed.
Using $v \cdot w = |v||w| \cos(\theta)$ and the fact that $a^2 + b^2 = r$, one should be able to set the two angles equal and solve for a single variable, say $a$.
I haven't explored just how harrowing the algebra becomes, but it certainly seems like a viable method.
A: 
Let the center of the circle be $O$.
Let $\phi = \angle AOB , \ \Psi = 90^\circ + \alpha = 90^\circ + \beta , \ \theta = \angle AOP $
Let $\overline{OA} = a , \ \overline{OB} = b , \overline{OP} = r $
Using the law of sines applied to $\triangle AOP $
$ \dfrac{a}{\sin \Psi} = \dfrac{r}{\sin( \theta + \Psi)} $
So that,
$ a ( \sin \theta \cos \Psi + \cos \theta \sin \Psi ) = r \sin \Psi $
from which
$ \tan \Psi = \dfrac{ a \sin \theta }{ r - a \cos \theta } $
Similarly, applying the law of sines to $\triangle BOP$ results in
$ \dfrac{b}{\sin \Psi} = \dfrac{ r }{\sin(\phi - \theta + \Psi ) } $
So that
$ b ( \sin(\phi - \theta) \cos \Psi + \cos(\phi - \theta) \sin \Psi ) = r \sin \Psi $
From which,
$ \tan \Psi = \dfrac{ b \sin(\phi - \theta) }{r - b \cos(\phi - \theta) } $
Hence, we now have
$\dfrac{ a \sin \theta }{ r - a \cos \theta } = \dfrac{ b \sin(\phi - \theta) }{r - b \cos(\phi - \theta) } $
from which
$a \sin \theta (r - b \cos(\phi - \theta)) = b (\sin(\phi - \theta) )(r - a \cos \theta) $
Multiplying this out:
$a r \sin \theta - a b \sin \theta \cos(\phi - \theta) - b r \sin (\phi - \theta) + a b \cos \theta \sin(\phi - \theta) = 0$
Expanding $\cos(\phi - \theta) $ and $ \sin(\phi - \theta)$ our equation becomes
$(- b r \sin \phi) \cos \theta + (a r + b r \cos \phi) \sin \theta + (a b \sin \phi) \cos 2 \theta + (- a b \cos \phi) \sin 2 \theta = 0 $
Which is of the form
$c_1 \cos \theta + c_2 \sin \theta + c_3 \cos 2 \theta + c_4 \sin 2 \theta + c_5 = 0 $
where $c_1, c_2, c_3, c_4, c_5$ are known constants.  This equation can be solved in closed form or numerically.  For the closed form solution check here
Once $\theta$ is determined, and to determine point $P$, find the cross product $\vec{OA} \times \vec{OB} $, with the $z$-component of both vectors set to zero.  If the cross product has a positive $z$-component, then $\vec{OB}$ is counter clockwise from $\vec{OA}$  otherwise it is clockwise.  Now $P$ is a rotation of $\bigg(r \dfrac{\vec{OA}}{a}\bigg)$ by a positive $\theta$ in the first case, and by a negative $\theta$ in the second case.  Rotation is achieved using the rotation matrix
$ R = \begin{bmatrix} \cos \theta && - \sin \theta \\ \sin \theta && \cos \theta \end{bmatrix} $
So,
$ P = \dfrac{r}{a} R A $
