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}$?

 A: 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).
A: This adds on to @kyary 's answer.
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.

*

*UPDATE: Here is @kyary 's equation [already in Snell Law form]
 solved by WolframAlpha (be patient):
https://www.wolframalpha.com/input/?i=0%3Dx%2Fsqrt%28d%5E2%2Bx%5E2%29-%28a-x%29%2Fsqrt%28c%5E2%2B%28a-x%29%5E2%29R


0=x/sqrt(d^2+x^2)-(a-x)/sqrt(c^2+(a-x)^2)R
 where $R=n_2/n_1$.

Without starting explicitly from Snell's Law...

*

*One can give the Fermat's Principle problem to WolframAlpha to solve symbolically (be patient):
https://www.wolframalpha.com/input/?i=0%3Ddiff%28+sqrt%28d%5E2%2Bx%5E2%29%2BRsqrt%28c%5E2%2B%28a-x%29%5E2%29%2Cx%29


0=diff( sqrt(d^2+x^2)+R*sqrt(c^2+(a-x)^2),x)
 where $R=n_2/n_1$.



*Here is an interactive visualization implemented in Desmos.
https://www.desmos.com/calculator/d21yvrq9io

For the given points with $n_1=1$ and $n_2=1.5$, we have the image below.
The visualization calculates $n_1\sin\theta_1=0.806200908137 =n_2\sin\theta_2$.


