Solving $P_y\frac{\sin{\theta_1}}{\sqrt{1-\sin^2{\theta_1}}}+F_y\frac{m\sin{\theta_1}}{\sqrt{1-m^2\sin^2{\theta_1}}}=F_x - P_x$ for $\sin\theta_1$

I am simulating the phenomenon of refraction, ruled by Snell's law: $$\sin\theta_1\cdot n_1 = \sin\theta_2\cdot n_2$$

For this question, is enough to know that relationship between both angles. $$n_1$$ and $$n_2$$ are given constants. The thick black line is called "interface".

What I am trying to achieve:

Based on user input, I know the coordinates of the point $$P$$ and point $$F$$ is fixed. Now, to draw the lines shown above, I need to determine where is the point $$I$$.

Specifically, I need to know its $$x$$ coordinate: the horizontal distance between $$P$$ and $$I$$.

How I tried to do it:

Since I know the horizontal distance between $$P$$ and $$F$$, I can equate it to the sum of the opposing legs of two right triangles, formed by angles $$\theta_1$$, $$\theta_2$$ and the vertical thin line, as shown below:

$$P_y\cdot\tan{\theta_1}+F_y\cdot\tan{\theta_2}=F_x - P_x$$

(Assuming absolute values and the interface is located at $$y=0$$).

Using Snell's law, I found that $$\sin{\theta_2}=\frac{n_1}{n_2}\sin{\theta_1}$$. Then, I tried to substitute it in the equation above. To simplify the expression, let $$m=\frac{n_1}{n_2}$$.

The issue

In terms of $$\sin{\theta_1}$$, the equation becomes:

$$P_y\frac{\sin{\theta_1}}{\sqrt{1-\sin^2{\theta_1}}}+F_y\frac{m\sin{\theta_1}}{\sqrt{1-m^2\sin^2{\theta_1}}}=F_x - P_x$$

And I'm stuck here. WolframAlpha could not help, and I don't know any numerical method good enough to determine $$\sin{\theta_1}$$. How can I solve that equation for $$\sin{\theta_1}$$?

• Maybe try midpoint bisection on $I_x$, bounded at first by $P_x$ and $F_x$, replacing one or other of the bounds if $n_1 \sin\theta_1 < n_2\sin \theta_2$ or not. May 3 '21 at 20:12
• Very good use of computational resources. By making an educated solution (say, always $0.75\cdot (P_x - F_x)$) it should converge pretty quickly. I probably will go with @kyary's solution though, even if it requires me to use the Newton–Raphson method, it will be computationally cheaper than many IFs for the bisection. Cheers! May 3 '21 at 20:30
• Determining the point I is essentially solved by Fermat's Principle of Least Time for refraction ( feynmanlectures.caltech.edu/I_26.html#Ch26-S3 ), which leads to the algebraic result in rectangular coordinates given by @kyary). May 3 '21 at 21:06
• It might be more useful if title includes "Refraction" or "Snell's Law" May 3 '21 at 22:48

There's an easier way to do it, assuming you have a way of solving quartics. By the way quartics can be solved systematically.

Plugging our values into snell's law, we get

$$\frac{x}{\sqrt{d^2 + x^2}}\cdot n_1=\frac{a-x}{\sqrt{c^2+(a-x)^2}}\cdot n_2$$

Squaring both sides gives us a quartic. Assuming you solve for x, you can of course use it to find sin(theta 1).

As I mentioned in the comments to the OP, one could use Fermat's principle of least time for refraction. (See Fig 26.4: https://www.feynmanlectures.caltech.edu/I_26.html#Ch26-S3 )

The optical path length is [using the variables given by @kyary] $$d_{opt}=n_1\sqrt{d^2+x^2}+n_2\sqrt{c^2+(a-x)^2}.$$ One should find the critical value of $$x$$ that extremizes this.

The result is Snell's Law, as expressed in @kyary 's answer.

Without starting explicitly from Snell's Law...