Reflect a ray off a circle so it hits another point

Solution 1:

Made a program which seems to work very well, i.e. it has passed all the examples I have tried. It uses no iteration but does require solving a cubic polynomial. I'll provide the code if you are interested.

My idea is the following: If I rotate the points E and P around the circle center, then at some rotation angle $\alpha$ the point of reflection must be at the top of the circle, i.e. at the point $(0,1)$. And a way to know I have found the right rotation angle is if the slope of a line from E' to $(0,1)$ is the negative equal of the slope from P' to $(0,1)$. If I could find the angle in this way, then the point of reflection, in the original setup, would just be located at an angle of $\frac{\pi}{2} - \alpha$.

So, if we let $E=(a,b)$ and $P=(c,d)$, then the rotated position of $E$, call it $E'=(a_1, b_1)$, would be:$$a_1 = a \cos(\alpha)-b \sin(\alpha)$$ $$b_1 = a \sin(\alpha)+b \cos(\alpha)$$ Same goes for $P'=(c_1,d_1)$.

If the slopes from $E'$ and $P'$ to (0,1) are negative equals, then $$\frac{b_1-1}{a_1} = -\frac{d_1-1}{c_1}$$

Inserting the equations for $a_1, b_1,$ etc into the equation above and simplifying, we get $$(ac-bd)\sin(2\alpha) + (bc+ad) \cos(2\alpha)+(b+d) \sin(\alpha) -(a+c)\cos(\alpha) = 0$$ Not so easy to solve. But, I decided to try an approximation, simply substituting the trig functions with the first two terms of their Taylor series. These are quite accurate up to about $40^\circ$, which is enough as I'll later talk about.

Inserting the substitutions and simplifying we get $$-\frac{8e+g}{6}\alpha^3 + (\frac{h}{2}-2f)\alpha^2+(2e+g)\alpha+(f-h) = 0$$ where $e=(ac-bd), f=(bc+ad), g=(b+d)$ and $h=(a+c)$.

The roots of this cubic equation can now be found in the usual manner.

This was my first version of the program and it worked somewhat, but could be off by $5-6 ^\circ$ which is too much. The problem is that the approximations get inaccurate for large angles. The solution I found was to do an initial rotation of E and P such that the angle bisector of $\angle E'CP'$ is at $90^\circ$. This puts E and P close to where they should be and means the final adjustment angle to be found via the cubic equation, is small.

Based on the examples I've tried (and checked via GeoGebra drawings), the accuracy of the found angle is now close to 2 decimal places. An example below:

The input was $E=(-1,-1.5), P=(-2,2)$ and the ouput was the angle $-168.05^\circ$.