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 Schrödinger 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 3 P0,2 states as well as the 1 S0. The two-gamma quarkonium decays that we investigate are for the η c , η ′ c , χ c0 , χ c2 , π 0 , π2, a2,and f ′ 2 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 π 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.