Express generic "penthagorean" triplet If asked to express the generic pythagorean triplets satisfying 
$a^2+b^2=c^2$ you would answer 
$$
a = k(r^2 - s^2), b = 2krs, c = k(r^2+s^2)
$$
with $k,r,s \in \Bbb N$ and $r>s$  and $r \ne s \mod 2$ and $\gcd(r,s) = 1$.  This generates every pythagorean triplet exactly once as the pair $(r,s)$ takes on its allowed values.
I'm looking to derive an analogous formula to generate every "penthagorean" triplet, that is, every instance where the sum of two pentagonal numbers is a pentagonal number.  That is, find all $(a,b,c)$ satisfying
$$
\frac{3a^2-a}2+\frac{3b^2-b}2=\frac{3c^2-c}2
$$
Deriving the formula for pythagorean numbers, using only elementary number theory such as considerations in mods $2$ and $4$ is so straightforward that I thought the analogous formula would be easy to obtain, but I am getting stuck.
 A: Here are a few one-dimensional parametric forms, if it helps. This was found by assuming the solution space was a line in $\mathbb{N}^3$.
$$
\begin{align}
&(16x+4, 30x+7, 34x+8) \\
&(33x+4, 56x+7,65x+8) \\
&(7x+7, 24x+23, 25x+24) \\
\end{align}
$$
A: from pages 124,125 of Magnus. If we have integers $a,b,c,d$ with $ad-bc=1$ and $a+b+c+d \equiv 0 \pmod 2,$ and we take
$$
\left(
\begin{array}{ccc}
\frac{1}{2} \left( a^2 + b^2 + c^2 + d^2 \right)&ab+cd&\frac{1}{2} \left( a^2 - b^2 + c^2 - d^2 \right) \\
ac+bd&ad+bc&ac-bd \\
\frac{1}{2} \left( a^2 + b^2 - c^2 - d^2 \right)&ab-cd&\frac{1}{2} \left( a^2 - b^2 - c^2 + d^2 \right) \\
\end{array}
\right)
\left(
\begin{array}{c}
x \\
 y\\
z \\
\end{array}
\right) =
\left(
\begin{array}{c}
u \\
v \\
w \\
\end{array}
\right)
$$
THEN
$$ u^2 - v^2 - w^2 = x^2 - y^2 - z^2  $$
At a minimum, the modular group takes a solution to your problem to another. Oh, multiply by 12 and complete the square,
$$  (6p-1)^2 + (6q-1)^2 = 1 + (6r-1)^2  $$
so that $6r-1$ becomes either $x$ or $u$
The observations in Magnus (1974) go back to Fricke and Klein (1897)
I cannot recall whether this gives all solutions to $x^2 - y^2 - z^2 = -1.$ Worth experimenting. The 3 by 3 matrix has determinant 1 and should preserve $\gcd(x,y,z).$
A: A $3$-parameter formula which always generates triplets of the form $$
x^2+y^2=z^2+1
$$ 
is (with $m,n,p \in \Bbb{Z}$)

$$
x = 2(mn+1)(mp-1)+1 \\
y = pm^2(n^2-1)+m(2np-n^2+1)+p-n \\
z= pm^2(n^2+1)+m(2np-n^2-1)+p-n
$$ 

I'm not certain that this formula generates all $(x,y,z)$ triples satisfying $x^2+y^2=z^2+1$ or even that it generates all such triples that can be formed using the Fricke/Klein technique, but numerical exploration indicates that it seems to cover them all.  Unlike the 3-parameter pythagorean formula, this formula does repeat triples, some in trivial ways but others by non-trivial coincidence.

Later observation - This 3-parameter formula does not generate all the triples satisfying $x^2+y^2=z^2+1$.

As noted by @Will Jagy, when all three elements $(x,y,z)$ of the triplet are of the form $6k-1$ you can form a penthagorean triplet $(\frac{r(3r-1)}2,   \frac{s(3s-1)}2, \frac{t(3t-1)}2)$ using
$$
(r = \frac{x+1}6, s =\frac{y+1}6, t = \frac{z+1}6) 
$$
