What we know: length of $a$, $b$, $c$
$a$ and $b$ are the sides of the large rectangle
what I need to find:
length of $d$ or just one angle (yellow color)
I found some solution on the internet with a similar problem, but damn is so hard to understand for me http://people.missouristate.edu/lesreid/POW10.html
I believe that is the way to calculate "4 collision places" in a more simple way
Let $\alpha$ be the angle with interrogation mark (that is found back in each of the NE, NW, SE, SW similar right triangles), you have the 2 following relationships between quantities $a,b,c,d,\alpha$:
$$\begin{cases} d \cos \alpha + c \sin \alpha &=& a \\ c \cos \alpha + d \sin \alpha &=& b\end{cases} \tag{1}$$
Said otherwise, and using the idea given by IAmAGuest in his comments (thanks to him), (1) can be considered as a linear system
$$\begin{cases}d x + c y &=& a \\ c x + d y &=& b\end{cases} \ \text{with} \ x:=\cos \alpha, y:=\sin \alpha\tag{1a}$$
The solutions of system (1a) are expressible in the following way (Cramer formulas):
$$x=\dfrac{\begin{vmatrix}a&d\\b&c\end{vmatrix}}{\begin{vmatrix}c&d\\d&c\end{vmatrix}}=\dfrac{ac-bd}{c^2-d^2}, \ \ \ y=\dfrac{\begin{vmatrix}c&a\\d&b\end{vmatrix}}{\begin{vmatrix}c&d\\d&c\end{vmatrix}}=\dfrac{cb-ad}{c^2-d^2},\tag{1b}$$
Constraint $x^2+y^2=1$ gives rise to the following relationship:
$$(ac-bd)^2+(cb-da)^2=(c^2-d^2)^2\tag{1c}$$
As a consequence, $d$ is "a"/"the" real positive solution of the 4th degree equation in variable $d$:
$$d^4-(a^2+b^2+2c^2)d^2+4(abc)d-(a^2+b^2-c^2)c^2=0\tag{1d}$$
Taking the numerical values $a=4,b=3,c=1$ considered in the (complicated) reference you gave, one obtains the same real positive root as theirs : $d\approx 4.0851$.
Remark: An alternative (longer) method is worth mentionning. It amounts to transform (1) into polynomial expressions using the classical rationalization formulas:
$$\begin{cases}\cos \alpha &=& \frac{1-t^2}{1+t^2} \\ \sin \alpha &=& \frac{2t}{1+t^2}\end{cases} \ \ \ \text{where} \ \ \ t:=\tan \frac{\alpha}{2}\tag{2}$$
Plugging expressions (2) into (1) generates 2 quadratic equations in variable $t$:
$$\begin{cases}2ct + d(1-t^2) &=& a(1+t^2) \\ c(1-t^2)+ 2dt&=& b(1+t^2) \end{cases} \ \iff \ \begin{cases}(a+d)t^2-2ct +(a-d)&=&0\\ (b+c)t^2-2dt+(b-c)&=& 0 \end{cases}\tag{3}$$
that must share a common root. This can be expressed by eliminating the $t^2$ term: indeed multiplying (3a) by $(b+c)$ and (3b) by $(a+d)$ and subtracting allows to express $t$ as a function of $a,b,c,d$. Plugging this expression into (3a) for example gives back (1d) .
If you are acustomed to "resultants" (see for example here), expressing that the two polynomials in (3) share a common root can be done by setting the following resultant to zero:
$$\begin{vmatrix}a+d&-2c&a-d&0\\0&a+d&-2c&a-d\\b+c&-2d&b-c&0\\0&b+c&-2d&b-c\\\end{vmatrix}=0. \tag{4}$$
Expanding (4) gives back (1d), once more.
Edit 1: I just found [using the formula searching engine approach0.xyz] this reference with the same numerical approximation.
Edit 2:: Applying Descartes' rule of signs to (1d) where the coefficients between parentheses are all $>0$, one can say that equation (1d) has at least $1$ positive root.