All the primitive Pythagorean Quadruples are known. This is Theorem 3 on page 176 and Theorem 4 on page 177 of Jones_Pall_1939.pdf, available at TERNARY as a pdf. The same information is on the first two pages of Pall_Automorphs_1940.pdf at the same site.
The short version is this: you have an odd number $ W $ Find all quadruples $a,b,c,d$ with $$ a^2 + b^2 + c^2 + d^2 = W, $$ where we are allowed to mix order, take the variables to be positive, negative, or zero. Then all the primitive quadruples, odd entry first, are given by
$$ \left( a^2 + b^2 - c^2 - d^2 \right)^2 + 4 \left( ad-bc \right)^2 + 4 \left(ac+bd \right)^2 = W^2. $$
Meanwhile, if $t$ is some (positive) divisor of $W,$ we can do the same process for $W/t$ and $W^2 / t^2,$ then multiply all $a,b,c,d$ by $t.$
July 2015: for another project, I decided to try to generate all the quadruples by the sort of three-parameter formulas one gets by stereographic projection to $\mathbb S^2.$ the results are really disappointing. I think I will stick with the four parameter thing. I found it in Jones and Pall, but it goes back at least to V. A. Lebesgue (and likely known to Euler), https://en.wikipedia.org/wiki/Pythagorean_quadruple#Parametrization_of_primitive_quadruples
and is a simple calculation using quaternions with ordinary integer coefficients. Alright, the (first correct) proof that the four parameter recipe gives all primitive solutions is from 1920, by Dickson. In 1941, Skolem gave a proof that can be adjusted to give a reasonably direct algorithm for taking a quadruple and reconstructing the four parameters. In 1956, one F. Steiger gave inequalities that make the mapping one-to-one. This is all reported in a 1962 article by Robert Spira, in the maa Monthly, May 1962, volume 69, number 5, pages 360-365, title The Diophantine Equation $x^2 + y^2 + z^2 = m^2.$ One thing I did not initially notice, Spira just discards the Pythagorean triples; in the quadruple setting, if one of the numbers is $0,$ we have a triple and can easily recover parameters.
Let's see, I was able to adapt the methods of Jones and Pall to $x^2 + y^2 + z^2 = 3 t^2,$ coefficients other than $3$ could also be handled this way. http://math.stackexchange.com/questions/1964607/when-will-a-parametric-solution-generate-all-possible-solutions/1965805#1965805
Every square can be written as the difference of two squares.
$(x+1)^2-x^2=2x+1$
This shows that any odd number, in particular an odd square, is the difference of two squares.
$(x+1)^2-(x-1)^2=4x$.
Any multiple of 4 is the difference of squares. An even square will be divisible by 4, so even squares are also the difference of two squares.
Best Answer
WLOG for $d^2=a^2+b^2, a=2mn, b=m^2-n^2$
$a+b=m^2-n^2+2mn$
Let $m^2-n^2+2mn=(m+pn)^2\iff n(p^2+1)=2m(1-p)$
$\dfrac n{2(1-p)}=\dfrac m{1+p^2}=r$(say)
Now $(1+p^2,1-p)|(1+p^2,1-p^2)$ which must divide $1+p^2+1-p^2=2$
Case $\#1:$ $(1+p^2,2(1-p))=1$ if $p$ is even
$n=2r(1-p^2),m=r(1+p^2)$
Case $\#2:$ For odd $p, p^2\equiv1\pmod8, 1+p^2\equiv2$
In that case, $(1+p^2,2(1-p))=2$
$n=r(1-p)$ and $m=\dfrac{r(1+p^2)}2$