Solutions to the Diophantine equation $x^2+3y^2+3z^2=n$ For a fixed positive integer $n$, the Diophantine equation
$$x^2 + y^2 + z^2 = n$$
was studied by Gauss in Disquisitiones Arithmeticae. As is known, this equation is intimately connected to the quaternion algebra $B_{-1,-1}$, characterized by the relations
$$
i^2 = -1, j^2 = -1, k^2 = (ij)^2 = -1.
$$
More precisely, it is equivalent to the norm equation
$$
\operatorname{nrd}(xi + yj + zk) = n,
$$
where $\operatorname{nrd}$ denotes the reduced norm of the (integral) pure quaternion $xi + yj + zk$. For me, the most interesting fact here is that for $n \not \equiv 0 \pmod 4$ the number of solutions to this equation is some multiple the Hurwitz class number $H(n)$. In the special case $n \equiv 7 \pmod 8$, it is a multiple of zero.
Now, I wonder if the same phenomenon was observed for other Diophantine equations of this kind. For example, in the quaternion algebra $B_{-1, -3}$ defined by the relations
$$
i^2 = -1, j^2 = -3, k^2 = (ij)^2 = -3,
$$
we can consider the norm equation
$$
\operatorname{nrd}\left(ix + \frac{1+j}{2}y + \frac{i + k}{2}z\right) = n,
$$
which is equivalent to the Diophantine equation
$$
x^2 + y^2 + z^2 + xz = n.
$$
Is there a formula for the number of its solutions? (see the OEIS sequence A014453) Is it true that, for some $n$, it is equal to the multiple of the (Hurwitz?) class number of some number field, say $\mathbb Q(\sqrt{-n})$? Perhaps, a Diophantine equation
$$
x^2 + 3y^2 + 3z^2 = n
$$
was studied in detail? I'd be thankful for any references.
 A: One of the beautiful (and sometimes flummoxing) aspects of the theory of ternary quadratic forms is that you can find answers and questions coming from many different points of view!  Using modular forms will give you an answer along the lines that Henri Cohen suggests.  I'd like to put in a pitch for a quaternionic approach, which has other features.
A reference is my book http://quatalg.org.  The theorem of Gauss is Theorem 30.1.3.  You can mimic the proof of this theorem for your ternary quadratic forms, to express the number of primitive representations in terms of (Hurwitz) class numbers, because the associated quaternion orders have type number one.  
Here are a few more details--if you'd like a more complete write up, just ask.  The main input is Theorem 30.4.7, which sums the number of optimal embeddings = number of primitive representations over the right class set of an order, expressing this in terms of the class number and a local correction factor.  The two quadratic forms you list arise as the reduced norm on the trace zero submodule of orders $\mathcal{O}$ in the quaternion algebra $B=(-1,-3\,|\mathbb{Q})$ of discriminant $3$.  The first form $x^2+y^2+z^2+xz$ arises from $\mathcal{O}$ a maximal order one: this order is in fact Euclidean under the reduced norm, and so has class number $1$, which means the average is just over the one term.  The second form $x^2+3y^2+3z^2$ arises from an order of reduced discriminant $36$ and index $12$; I compute using Magma that this order has class number $2$ but type number $1$, which means that the sum is over two terms that are equal, so we just need to divide the right-hand side by $2$.  
So to finish, we need to compute the local embedding numbers, which provides a correction term depending on $n$ modulo a power of $2$ and $3$.  This gets a bit technical, so I would just stare at the numerical answer provided by Will Jagy.  For the first form, the order is maximal so the local embedding numbers are in 30.5 (just a correction at $3$); unfortunately for the second form requires a separate calculation (not covered by 30.5 or 30.6).
Happy to say more if you'd like!
A: The theta function associated to your quadratic form (here $\theta(\tau)\theta(3\tau)^2$) belongs to the modular form space $M_{3/2}(\Gamma_0(12))$,
and you are in luck because there are no cusp forms, so the space is entirely
generated by three Eisenstein series whose coefficients are similar to Hurwitz
class numbers and easily computable, so there is an explicit formula analogous
to Gauss in your case (which I did not take time to work out). This works
as long as there are no cusp forms. However, for instance $S_{3/2}(\Gamma_0(28))$ is nonzero, so you could not do the same with $x^2+7y^2+7z^2$.
