I believe I've found a solution that is quite basic and should be approachable by people with an elementary understanding of mathematics. Color each square of the grid with colors $1$ through $k$, as follows. In the image below, $k=7$.
$$\begin{array}{|c|c|c|c|c|c|c|c|c|c|c}
\hline
1&2&3&4&5&6&7&1&2&3&\dots\\\hline
2&3&4&5&6&7&1&2&3&4&\dots\\\hline
3&4&5&6&7&1&2&3&4&5&\dots\\\hline
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\
\end{array}$$
Then each tile will contain exactly one color of each kind. It follows that
If the grid can be tiled by $k\times 1$ rectangles, then it contains the same number of squares of each color.
Let's agree that the grid is $m$ rows by $n$ columns.
If $m$ is divisible by $k$, say $m=ck$, then for each color and each column there will be exactly $c$ squares of that color in that column. Hence, the total number of squares of each color is the same. Something similar holds when $n$ is divisible by $k$.
Now, suppose neither $m$ nor $n$ is divisible by $k$, and write $m=ck+r$ and $n=dk+q$, with $1\leq r,q\leq k-1$. The subgrid $ck\times n$ contains an equal number of squares of each color, by the previous argument. Consider then the bottommost strip, of size $r\times n$.
Once again, the subgrid of size $r\times dk$ contains an equal number of squares of each color. We are finally left with a subgrid of size $r\times q$.
Assume without loss of generality that $r\leq q$. Renumbering if needed, it should look something like this:
$$\begin{array}{|c|c|c|c|}
\hline
1&2&\dots&q\\\hline
2&3&\dots&q+1\\\hline
\vdots&\vdots&\ddots&\vdots\\\hline
r&r+1&\dots&r+q-1\\\hline
\end{array}$$
Since $r\leq q$, it's easy to see that color $q$ shows up in exactly $r$ squares -- once on each row (consider the antidiagonal). On the other hand, color $1$ definitely does not show up on the second row, so it can appear at most $r-1$ times.
This shows that, as a whole, the $m\times n$ grid does not contain the same number of squares of each color. Therefore, if $k$ divides neither $m$ nor $n$, the $m\times n$ grid cannot be tiled by $k\times 1$ rectangles.
The problem in your code is
var nextVertex = intersectionPt
which should really be
var nextVertex = [0, 0]
(or at least a constant for all the tiles you generate).
The formula $\sum_{j=0}^4 n_j \mathbf e_j$ for the corner locations produces coordinates in a fixed coordinate system -- it's not relative to the intersection you start out with!
So the gaps you're seeing correspond to the distances between intersection points in the pentagrid, which you're erroneously adding to the tile coordinates.
A different problem that might bother you after you fix that is that the tiles you generate are too large compared to your pentagrid. You seem to be making the spacing between neighboring pentagrind lines equal to side length of the rhombs -- but the average distance between the "ribbons" corresponding to pentagrid lines is clearly longer than that. The tile side length is the width of a ribbon, and there's plenty of space left over between the ribbons!
(Image from https://www.ams.org/publicoutreach/feature-column/fcarc-ribbons)
This is not a problem for generating an infinite tiling, but it will make it difficult to predict which part of the pentagrid to look for intersections in, if you want to find all the tiles in a particular area of the coordinate system.
Curiously, none of your references seem to mention this difference in scale explicitly -- those that give actual formulas seem to assume that the lines on the pentagrid are unit spaced and that the rhomb sides are a unit apart. However, those references also don't actually show the complete pentagrid in the same plane as the finished Penrose tiling. Apparently in their view the pentagrid is only there to define which coordinate sets in $\mathbb Z^5$ correspond to corner points in the tiling, but the actual tiles live in a completely different coordinate system.
If you want each tile corner to appear roughly in the vicinity of the region of the pentagrid that corresponds to it, you need to scale one of the coordinate systems to match. Your sources don't give you the scaling factor to use, so let's derive it. We'll keep the tile side lengths as $1$, and figure out how far apart the pentagrid lines must be.
Suppose we're in the region of the pentagrid in your question with 5-coordinate tuple $(0,0,0,0,0)$, and we now move some large distance $D$ due right. What is the five-tuple corresponding to the region we land in, if the grid spacing is $G$?
Well, we move past $\frac DG$ red lines, crossing them at right angles. The orange lines make an angle of $72^\circ$ with the red ones, so the number of them we cross is an integer close to $\frac DG \cos 72^\circ$. The rounding becomes negligible when we let $D$ got to infinity, so just pretend it is $\frac DG\cos 72^\circ$ exactly. Similarly, the green lines make an angle of $144^\circ$ to the red ones, so the number of them we cross is $\frac DG \cos 144^\circ$. (The factor $\cos 144^\circ$ is negative, corresponding to the fact that we cross these lines in the direction of decreasing region coordinates). And so forth with the cyan and blue, with angles of $216^\circ$ and $288^\circ$.
All in all, we end up with the 5-coordinates, up to rounding
$$ (n_0,n_1,n_2,n_3,n_4) \approx \bigl(\tfrac DG,\tfrac DG\cos 72^\circ,\tfrac DG\cos 144^\circ,\tfrac DG\cos 216^\circ,\tfrac DG\cos 288^\circ\bigr). $$
The coordinates of the tile corner corresponding to that region is
$$ \sum_{j=0}^4 n_j \mathbf e_j = \sum_{j=0}^4 n_j \bigl(\cos(j\cdot 72^\circ), \sin(j\cdot 72^\circ)\bigr). $$
We're just interested in the $x$-coordinate, so plug in the $n_j$s above to get
$$ x \approx \frac DG \sum_{j=0}^4 \cos^2(j\cdot 72^\circ)
= \frac DG \sum_{j=0}^4 \frac{1 - \cos(j\cdot 144^\circ)}2
= \frac DG \cdot \frac{5-\sum_{j=0}^4 \cos(j\cdot 144^\circ)}2
= \frac DG \cdot \frac52 .
$$
(The last equals sign is because the sum of cosines is the change in $x$-coordinate from going all the way around a pentagram with unit sides).
If this $x$ is to grow at the same average rate as $D$, we must set
$$ G = \frac52. $$
Beware, however, that even if you scale the grid by this factor, you still won't get closer than the tiles being "roughly in the vicinity" of the grid intersections that represent them. The nice images in the presentations you link to make it look like each rhombus will contain the grid intersection it comes from -- but it's not possible to align the grid such that this is always the case, nor such that every tile corner is inside the grid region that you compute its location from.
The pentagrid contains places where five grid lines almost but not quite meet in a point, like in this region of your pentagrid:
In this particular almost-intersection the fit is loose enough that we can see it's not exact -- but if your continue the pentagrid far enough, you will eventually find places where almost-intersections like this take place within arbitrarily small areas. The configuration has 6 internal regions and 10 actual intersections, and each of those intersections generates a full-sized tile in the final tiling. The 10 tiles come together in this pattern:
whose diameter is somewhat larger than the pentagrid spacing! Clearly not all of the tiles can lie directly above their defining intersection.
(We can see a differently oriented instance of this pattern in your buggy picture. Those tiles have particularly small gaps between them because the intersections that represent them are close together).
Best Answer
A recent paper by Lutfalla and Fernique provided one possible answer to this question. The restrictions imposed by the possible vertex figures do not suffice, as demonstrated in this counterexample:
But if we require the tiles touching every edge emanating from a given vertex to come from a valid Penrose tiling, i.e. one of the following 15 configurations:
then it is forced to be a valid Penrose tiling.
This is not that compactly describable, though, so I'm not accepting this self-answer yet in hopes that a better classification of unmarked Penrose tilings is posted.
Edit 2023-05-14: Theorem 6.1 of M. Seneschal's Quasicrystals and Geometry states that any tiling whose vertex figures ("stars") are one of the seven allowable configurations and such that no two stars that share a rhomb are related by a half-turn about the center of that rhomb is a Penrose tiling. (Note that this second condition fails to obtain in the above counterexample - there are two kinds of vertices with the aforementioned 180 degree symmetry.)