This is just elementary comment, but may be too long for comment, probably you know it, but just for completeness.

I happened to ask almost (but not exactly) the same question some days ago: 

https://mathoverflow.net/questions/104439/conjugcy-classes-in-glf-2-glf-q

It sounds a little different - size of conjugacy classes in GL(F_q).
But conjugacy class of element "C" is size of GL(F_q)/ size of centralizer of element "C".
Because size of any orbit is |G|/|Stabilizer| and here we have action by conjugation and stabilizer
is centralizer.
So if you know the size of centralizer - you know size of conjugacy class.

The difference is, of course, that you asked about Mat(F_q) while setup above is about non-degenerate matrices GL(F_q).

-------------

Let me also write down some elementary facts for completeness.

If you consider element $C$ such that its characteristic polynomial is irreducible, (then it automatically minimal), then size of centralizer in Mat_n(F_q) is q^n (if I understand correctly)
and q^n-1 in GL_n(F_q).

One has good way to think about this centralizer: let consiser p(x) - char.pol. of "C".
F_q[x]/p(x) is a field F_q^n.
To any element of the field "a" in F_q^n one can correspond a matrix M(a) - the matrix of multiplication by "a":F_q^n -> F_q^n.

Map "a->M(a)" is clearly homomorphism of algebras so in particular all matrices M(a) commute among themselves. They provide a centralizer of a matrix M(x).