  Two-Body Dirac equations of constraint dynamics provide a covariant framework
to investigate the problem of highly relativistic quarks in meson bound states.
This formalism eliminates automatically the problems of relative time and
energy, leading to a covariant three dimensional formalism with the same number
of degrees of freedom as appears in the corresponding nonrelativistic problem.
It provides bound state wave equations with the simplicity of the
nonrelativistic Schroedinger equation. Here we begin important tests of the
relativistic sixteen component wave function solutions obtained in a recent
work on meson spectroscopy, extending a method developed previously for
positronium decay into two photons. Preliminary to this we examine the
positronium decay in the 3P_{0,2} states as well as the 1S_0. The two-gamma
quarkonium decays that we investigate are for the \eta_{c}, \eta_{c}^{\prime},
\chi_{c0}, \chi_{c2}, \pi^{0}, \pi_{2}, a_{2}, and f_{2}^{\prime} mesons. Our
results for the four charmonium states compare well with those from other quark
models and show the particular importance of including all components of the
wave function as well as strong and CM energy dependent potential effects on
the norm and amplitude. The results for the \pi^{0}, although off the
experimental rate by 15%, is much closer than the usual expectations from a
potential model. We conclude that the Two-Body Dirac equations lead to wave
functions which provide good descriptions of the two-gamma decay amplitude and
can be used with some confidence for other purposes.
