How to determinate the 4×4 S matrix so that the P gets projected into Q, on the XZ (Y=0) plane?
Q = S P
[Math] Projection matrix to project a point in a plane
geometrymatricesprojective-geometry
Related Solutions
Much theory up front, the actual bug fixing is way down there. The theory covers the general approach of how you come up with such a projection matrix and servers as a general answer to this kind of question, while the bug fixing is pretty local and specific to your implementation.
Theory
The line
Suppose you have the (variable) coordinates of your input point $P$, and the (fixed) coordinates of your light source $L$. Then any point on the line between these two points can be expressed as a linear combination $\lambda P + \mu L$. And conversely any linear combination will represent a point on that line. So for the projected point $P'$ you get the condition $$P'=\lambda P+\mu L$$
The plane
Let $F\cdot Q$ represent the dot product between the coordinates of your plane $F$ and the coordinates of a point $Q$. A point $Q$ lies on that plane if and only if $F\cdot Q=0$. So for the projected point $P'$ you get the condition $$F\cdot P'=0$$
The intersection
Now you have to combine these two: find coefficients $\lambda$ and $\mu$ such that the second equation holds. There is a very elegant way to do this, known (to some) as Plücker's $\mu$.
\begin{align*} \lambda &= F\cdot L \\ \mu &= -F\cdot P \\ P' &= (F\cdot L)P - (F\cdot P)L \end{align*}
The point chosen in this fashion will automatically lie on the plane, therefore will be your desired projected point.
As a matrix
Now you simply have to split that definition into its individual components.
\begin{align*} P'_i &= \sum_{j=1}^4 F_jL_jP_i - \sum_{j=1}^4 F_jP_jL_i & \forall i&\in\{1,2,3,4\} \\ M_{ij} &= -F_jL_i + \begin{cases} F\cdot L &\text{if }i=j \\ 0 &\text{otherwise} \end{cases} & \forall i,j&\in\{1,2,3,4\} \\ M &= \begin{pmatrix} \lambda-F_1L_1 & -F_2L_1 & -F_3L_1 & -F_4L_1 \\ -F_1L_2 & \lambda-F_2L_2 & -F_3L_2 & -F_4L_2 \\ -F_1L_3 & -F_2L_3 & \lambda-F_3L_3 & -F_4L_3 \\ -F_1L_4 & -F_2L_4 & -F_3L_4 & \lambda-F_4L_4 \end{pmatrix} & \text{with }\lambda&=F\cdot L \end{align*}
Your code
So your matrix does look reasonable. Your dot corresponds to my $\lambda$, and all the rest looks sane as well. Your problem lies with the pointFmatrixF method. There you incorrectly assume that the last coordinate of the resulting point will be 1. In this case it will be 40, so after fixing that line you get projectedPoint: (40.0, 0.0, 0.0, 40.0) which is just another homogenous representantion of the point (1, 0, 0). Simply divide all other coordinates by the last coordinate.
Also notice that your pointFmatrixF does implicitely assume point[w] == 1. This is the case for your special input, but not true in general. So you probably should multiply by the fourth component as well.
Your question claims a light source at (0, 0, 10) which would be homogenized to (0, 0, 10, 1). Your code however has a light at (0, 0, 10, 0), which is a light infinitely far away in the z direction. Fixing that as well, you get projectedPoint: (40.0, 0.0, 0.0, 20.0) which is the (2, 0, 0) you expected.
I'll start by looking about how you can describe the projection of an arbitrary point $Q$ onto your plane, then derive a matrix notation from that.
Projecting a single point
The point $Q$ and its projection form a line which is perpendicular to the plane. All these perpendicular lines meet in a point $F=(A,B,C,0)^T$, which is the point at infinity orthogonal to your plane. A generic point on the line connecting $Q$ and $F$ can be described as a linear combination $Q+\lambda F$. (This inhomogenous formulation using a single parameter excludes $F$ itself, which is all right since $F$ will not be the projected point in any case.) Now you are looking for the point where this line connecting $F$ and $Q$ intersects the plane, i.e. $\left<Q+\lambda F,h\right>=0$, where $h=(A,B,C,D)^T$ is the vector describing your plane. This you can solve for $\lambda$:
$$\left<Q+\lambda F,h\right>=\left<Q,h\right> + \lambda\,\left<F,h\right>=0\\ \lambda=-\frac{\left<Q,h\right>}{\left<F,h\right>}$$
To avoid the division, you can also use a multiple of that point:
$$Q+\lambda F\;\sim\; \left<F,h\right>\,Q - \left<Q,h\right>F $$
This is the projection of a generic point $Q$.
Finding the matrix
If you assume the coordinates of $Q$ to be variables, then the coordinates of the result will be linear in these variables, so you can interpret the whole operation as a linear map and therefore write it as a matrix. You can also find the formula of that matrix like this:
$$ \left<F,h\right>\,Q - \left<Q,h\right>F = \left<F,h\right>\mathbb 1\cdot Q - Fh^T\cdot Q = \left(\left<F,h\right>\mathbb 1 - Fh^T\right)Q \\ = \begin{pmatrix} B^{2} + C^{2} & -A B & -A C & -A D \\ -A B & A^{2} + C^{2} & -B C & -B D \\ -A C & -B C & A^{2} + B^{2} & -C D \\ 0 & 0 & 0 & A^{2} + B^{2} + C^{2} \end{pmatrix}\cdot Q $$
If anyone reading this post should be interested in central instead of orthogonal projection, simply use that center of projection as $F$, and the above formula will help you compute the matrix for that as well.
Best Answer
I've done done math and I have the solution now. This is simple similar triangles.
You have a light point L and a P point you want to project on the ground. If you draw the lines you find two triangles with the proportions. The red side divided by the green side is equivalent to the purple side divided by the pink side + green side like
$\frac{P_x - L_y}{P_y - L_y} = \frac{Q_x - L_x}{L_y}$
If you solve for $Q_x$ you have $Q_x = \frac{L_y \cdot P_x - L_x \cdot P_y}{P_y - L_y}$
The drawing is for 2D but don't forget we're on 3D so we need to solve for $Q_z$
$Q_z = \frac{-L_z \cdot P_y + L_y \cdot P_z}{P_y - L_y}$
We now want to find a matrix $M$ that gives the following result when multiplied by any point $P$ (The vertical axis is y and wee need a point on the floor so $y= 0$)
$ M *\begin{bmatrix}P_x\\P_y\\P_z\\1\\\end{bmatrix} = \begin{bmatrix}Q_x\\0\\Q_z\\1\\\end{bmatrix}$
The question here is to use homogeneous coordinates. If you look at $Q_x$ you have a division. There's no way to divide anything inside a matrix multiplication. So you use the homogeneous coordinate to get the result we want.
So we need a matrix $M$ multipled by any point $P$ that gives us $\begin{bmatrix} L_y \cdot P_x - L_x \cdot P_y\\ 0\\ -L_z \cdot P_y + L_y \cdot P_z\\ L_y-P_y\\ \end{bmatrix}$
and then we divide every row by $Py-Ly$ to get $\begin{bmatrix} \frac{L_y \cdot P_x - L_x \cdot P_y}{P_y-L_y}\\ 0\\ \frac{-L_z \cdot P_y + L_y \cdot P_z}{P_y-L_y}\\ 1\\ \end{bmatrix} $
$M = \begin{bmatrix} L_y & -L_x & 0 & 0\\ 0 & 0 & 0 & 0\\ 0 & -L_z & L_y & 0\\ 0 & -1 & 0 & L_y \end{bmatrix}$
Ps: I'm sorry if I'm not explaining the right way or if I'm making any mistake but the result is there.