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:
Conjugcy classes in GL(F_2)? GL(F_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).