Quadrilateral: How to calculate equal length ratios on the edges from a given length ratio on the diagonal? I have a quadrilateral.

The diagonal F divides the quadrilateral in two triangles. For each triangle all the side lengths as well as all the angles are known.
Furthermore a point P on the diagonal F is given by the ratio u/F.
I am looking for the points O and Q i.e. the ratios a/E and c/G so that a/b=c/d.
It is: a+b=E, u+v=F, c+d=G
What is the most efficient way to calculate this? (I need to use it in a computer program).
 A: Let's call the point where F and E meet R, and the point where F and G meet S. The angles $\rho:=\sphericalangle ORP$ and $\sigma:=\sphericalangle QSP$ are known, the angle $\phi:=\sphericalangle RPO = \sphericalangle SPQ$ is the one quantity we want to calculate.
By the law of sine in $\triangle ROP$ we have $\frac{a}{u}=\frac{\sin(\phi)}{\sin(\phi+\rho)}$. Similiarly, we get from $\triangle PSQ: \frac{d}{v}=\frac{\sin(\phi)}{\sin(\phi+\sigma)}$.
The equation you want to reach is
$$\frac{a}{E-a}=\frac{a}{b}=\frac{c}{d}=\frac{G-d}{d}$$
If you plug in the formulas for $a,d$ from above, you get an equation that contains, beside the unkown $\phi$, only the known values $E,G,u,v,\rho,\sigma$.
It might be possible to actually solve the equation for $\phi$, but I doubt it. What you can do however, since you are using a computer anyway, is to solve the equation numerically. That is easy in this case as both sides of the equation, when considered as functions of $\phi$, behave monotonically in opposite ways! 
If $\phi$ is $0$, then $O=R$ and $Q=S$ and $a=d=0$, so $\frac{a}{b}=0,\frac{c}{d} =\infty$. If $\phi$ increases, then $a$ and $d$ increase, while $b$ and $c$ decrease. This means $\frac{a}{b}$ will increase, while $\frac{c}{d}$ decreases.
Now just use $\phi=1°, 2°, \ldots$ and calculate the corresponding values $\frac{a}{b},\frac{c}{d}$. In the beginning, you will have $\frac{a}{b}>\frac{c}{d}$, just as for $\phi=0°$. But at some point, this will change, so you have a $\phi_0$ where still $\frac{a}{b}>\frac{c}{d}$ but then a $\phi_1=\phi_0+1°$ where $\frac{a}{b}<\frac{c}{d}$. Then you know that the value of $\phi$ that you seek must be between $\phi_0$ and $\phi_1$.
Now you can iteratively calculate $\phi_2:=\frac{\phi_0+\phi_1}{2}$, and check if for $\phi_2$ you get $\frac{a}{b}>\frac{c}{d}$ or $\frac{a}{b}<\frac{c}{d}$. If the former, continue with a new $\phi_0:=\phi_2$ ($\phi_1$ stays the same). If the latter, set new $\phi_1:=\phi_2$ ($\phi_0$ stays the same in this case).
This allows you to calculate the $\phi$ you want to arbitrary precisions. Since you are calculating the corresponding $a,b,c,d$ values for $\phi_0,\phi_1$ anyway, you can stop doing the iteration if the values of $a$ and $c$ for $\phi_0$ differ by no more than some accuracy that is good enough for you from the corresponding values for $\phi_1$.
A: Here are two solutions both based on coordinates. The first one is the simplest.
Have a look at the figure below, where the given quadrilateral is $AFEC.$ 


*

*First solution
Let the coordinates of the given points be $A(x_A,y_A),C(x_C,y_C),E(x_E,y_E),(x_F,y_F),P(x_P,y_P).$
Let $X(x_X,y_X),Y(x_Y,y_Y)$ be the coordinates of the looked for points.
Line segments $AF$ and $CE$ can be described in a parametric way  resp. as : 
$$\tag{1.1}X \in \text{segment AF}\ \ \iff \ \exists \alpha, \ 0 \leq \alpha \leq 1, \ s.t. \ \begin{cases}  x_X=x_A+\alpha(x_F-x_A)\\
y_X=y_A+\alpha(y_F-y_A)
\end{cases}$$
$$\tag{1.2} Y \in \text{segment CE} \ \iff \ \exists \beta, \ 0 \leq \beta \leq 1, \ s.t. \
\begin{cases}x_Y=x_C+\beta(x_E-x_C)\\
y_Y=y_C+\beta(y_E-y_C) \end{cases}$$
For example, for segment $AF$, you can check that, if $\alpha=0$, $X$ is $A$, if  $\alpha=1$, $X$ is $F$, and if  $\alpha=1/2$, $X$ is the midpoint of segment $AF$.
A classical alignment condition (http://mathworld.wolfram.com/Collinear.html) is as follows:
$$X,P,Y \ \text{aligned}  \ \iff \ \det \begin{pmatrix}x_X & x_P & x_Y \\
y_X & x_P & y_Y\\ 1 & 1 & 1\end{pmatrix}=0   \ \iff \ \det \begin{pmatrix}x_A+a(x_F-x_A)& x_P & y_A+a(y_F-y_A) \\
y_A+a(y_F-y_A)  & x_P & y_A+a(y_F-y_A) \\ 1 & 1 & 1\end{pmatrix}=0$$
for the same value of $a$. 
The last equation is equivalent, by subtracting column 2 to the other columns, to
$$\det \begin{pmatrix}(x_A-x_P)+a(x_F-x_A)&  (y_A-y_P)+a(y_F-y_A) \\
(y_A-y_P)+a(y_F-y_A)  & (y_A-y_P)+a(y_F-y_A)\end{pmatrix}=0$$
The development of this determinant will give rise to a quadratic equation in variable $a$, whose unique positive solution $a_1$ gives $X$ and $Y$, by replacing $\alpha$ and $\beta$ in (1.1) and (1.2) by $a_1.$


*

*Second solution
We will here need $B$, the intersecting point of lines $CE$ and $AF$.
Let us define: 
$$\tag{2.0}\begin{cases}E=eB+(1-e)C\\
F=fB+(1-f)A\\
P=p E+(1-p)A
\end{cases} \ \ \ \text{and} \ \ \ \begin{cases}X=(1-a)A+aB\\Y=cB+(1-c)C.\end{cases}$$
Then it suffices to take $a$ given by formula
$$\tag{2.1}a=\dfrac{a_1-\sqrt{a_2}}{2e}$$ where 
$a_1=f + p e - p f + p e f$ and
$a_2= p^2 e^2 (f+1)^2- 2p^2 e f^2 - 2p^2 e f - 4p e^2 f + 2p e f^2 + 2 p e f +(1-p)^2f^2.$
and $c$ given by:
$$\tag{2.2}c=\dfrac{a - a p - e p + a e p}{a - p}.$$
(see Matlab program below).
Proof of these formulas: let us use barycentric coordinates with respect to triangle $ABC$ ; I take the convention to place between square brackets the 3 barycentric coordinates of a point. 
Formulas (2.0) can be written in terms of barycentric coordinates:
$$E=[0,e,1-e], \ \ \ \ F=[1-f,f,0],$$
$$P=p[0,e,1-e]+(1-p)[1,0,0]=[1-p,ep,p(1-e)]$$
$$X=[1-a,a,0], \ \ \ \ Y=[0,c,1-c]$$
$$X,P,Y \ \text{aligned}  \ \iff \det \begin{pmatrix}1-a & 1-p & 0 \\
a & ep & c\\ 0 & (1-e)p & 1-c\end{pmatrix}=0.$$
from which (2.2) can be deduced.
Now, let us justify formula (2.2). We have the following equivalences:
$$\dfrac{a}{b}=\dfrac{c}{d} \ \iff \ \dfrac{a}{f-a}=\dfrac{c}{e-c} \iff 
 \ ae=fc $$
This last condition, taking into account (2), gives rise to a quadratic equation in variable $a$ whose solution is given by formulas (1).
Matlab program:


 A=[0;0];B=[5;0];C=[2;2];
 e=0.6;E=e*B+(1-e)*C;
 f=0.8;F=f*B+(1-f)*A;
 p=0.7;P=(1-p)*A+p*E;
 a1=f + p*e - p*f + p*e*f;
 a2=p^2*e^2*(f+1)^2 - 2*p^2*e*f^2 - 2*p^2*e*f - 4*p*e^2*f+2*p*e*f^2 + 2*p*e*f + (1-p)^2*f^2;
 a=(a1-sqrt(a2))/(2*e);
 c=(a - a*p - p*e + a*p*e)/(a - p);
 X=(1-a)*A+a*B;
 Y=c*B+(1-c)*C;




