prime order panmagic squares  

basic ingredients of panmagic squares  
Latin square generating formula LS(a) factor: m  3 
Latin squares obtained by formula  
LS(a): L_{i,j} = (a j + i) % m ; i,j ε [0,..,m1]; a = 2 .. m  2  
"intermediate" or most basic panmagic squares amount: (m4)(m3)/2 
S = m * LS(a) + LS(b) + 1; a < m/2; b != a  
the a is limited to the first half of the range since using the other half merely introduces a horizontal panflip variant (which we want to avoid), with b = a the square becomes irregular (which we also don't want) There are (m3)/2 a's and (m4) b's satisfying the conditions. 

(m4)(m3)/2  
5 7 11 13 
1 6 28 45 

most essential panmagic squares amount: [(m4)(m3)/2](^{m1}_{4}) 
S = m * LS(a)[p,q,r,s] + LS(b) + 1; a < m/2; b != a; p < {q,r,s}  
The number p is located just right of the top left '0', q just below, r the last number of the first column and s the last number of the firts row. The thus positioned numbers regulate whether the square is in normalized position by imposing the conditions: p < q (transposition) p < s (horizontal pan flip) q < r (vertical pan flip) however once established that p < (q,r,s), the numbers q,r and s can be freely permuted merely introducing squares which need a vertical pan flip to be put in normalized position it's easily recognized the these squares are independent of the other squares, this fact and the first two (q < r not considered of course) conditions above compensate the random positioning of the 4 selected numbers [p,q,r,s] (ie 4!/(2^{2}/3!) the selection of these four numbers (from m1 possible numbers) can of course be done in (^{m1}_{4}) ways. In the most essential square both sets (q,r,s) and the remaining (m5) numbers are put in their natural ordering by the here imposed definition. 

[(m4)(m3)/2](^{m1}_{4})  
5 7 11 13 
1 90 5,880 22,275 

normalized positioned panmagic squares amount: [(m4)(m3)/2](^{m1}_{4}) 3!(m5)!(m1)! 
S = m * LS(a)[p,q,r,s]_{=P(q,r,s)=P([1..m1]\{p,q,r,s})} + LS(b)_{=P((m1)+1)} + 1; a < m/2; b != a; p < {q,r,s} 

All panmagic squares now can be found by using digit changing permutations on both the high and the low component latin squares. Note that the permutation denoted by =P([1..m1]\{p,q,r,s}) is a permutation on the remaining numbers all three mentioned permutations are independent of one another introducing the mentioned factors 3! (m5)! and (m1)! onto the total amount of panmagic squares obtainable from described method, which (I believe) are all different by construction. 

[(m4)(m3)/2](^{m1}_{4})3!(m5)!(m1)!  
5 7 11 13 
144 777,600 92,177,326,080,000 2,581,228,494,028,800,000 

panmagic squares amount: [(m4)(m3)/2](^{m1}_{4}) 3!(m5)!(m1)!m^{2} = [(m4)(m3)/2](m!)^{2}/4 
S = [m * LS(a)[p,q,r,s]_{=P(q,r,s)=P([1..m1]\{p,q,r,s})} + LS(b)_{=P((m1)+1)} + 1]_{@[x,y]}; a < m/2; b != a; p < {q,r,s}; x,y ε [0..m1] 

Panrelocating the topleft corner '1' gives the total amount of panmagic squares  
[(m4)(m3)/2](^{m1}_{4})3!(m5)!(m1)!m^{2}  
5 7 11 13 
3,600 38,102,400 11,153,456,455,680,000 436,227,615,490,867,000,000 

panmagic squares [(m4)(m3)/2](m!)^{2} 
S = [m * LS(a)_{=P(m)} + LS(b)_{=P(m)} + 1];a < m/2; b != a;  
Allowing all digits to get changed on both components of course also gives all panmagic squares, the factor 4 discrepancy between this number and the previous is due to the fact that also reflectional variants are obtained by this method 

Selfcomplementairy panmagic squares [(m4)(m3)/2]((m1)!!)^{2}/4 
S = [m * LS(a)_{=P(a;m)} + LS(b)_{=P(b;m)} + 1];a < m/2; b != a; with P(a;m) permutations such that P[(a+1)i] + P[(a+1)(m1i)] = m 

In order for the complete square to fullfill the selfcomplementairy condition, ie S[i,j] + S[m1i,m1j] = m^{2}+1 it is neccessairy for each component to be selfcomplementairy. A little figuring showed me that this is established with digit changing permutations which is merely some permutation of a symmetric permutation given exactly by the above stated condition. For a prime order m the amount of symmetric permutations is given by (m1)!!, the deviding factor of 4 takes care of the reflectional variants the method introduces 

[(m4)(m3)/2](m1)!!/4  
5 7 11 13 
16 3,456 103,219,200 23,887,872,000 
prime order panmagic squares  

basic ingredients of panmagic squares  
pandiagonal latin square permutation generated LS(P[0..m1]) factor: F_{m} currently known: F_{13} = 348 
Latin squares obtained using the natural permutatio on the xaxis and a permutation on the yaxis satisfying the pandiagonal condition 

LS(P[0..m1]): L_{i,j} = (P[j] + i) % m ; i,j ε [0,..,m1] {P[i] + m +/i; i = 0..m1} ε {P[0..m1]} 

let P[0..m1] be a permutation of the numbers 0..m1 and {P[0..m1]} the set of all these permutations (m! members), with P_{m,0} (the natural permutation of m digits is [0,1,...,m1,m1] on the xaxis, and some other permutation on the yaxis, in order to obtain a pure latin square the diagonal and subdiagonal need to consitst of permutations themselves, this results in the two conditions: {P[i] + m +/i; i = 0..m1} ε {P[0..m1]} This ensures that the latin square: LS(P[0..m1]): L_{i,j} = (P[j] + i) % m ; i,j ε [0,..,m1] is a {pandiagonal} latin square 

Orthogonal pairing factor O_{m} yet to be determinined 
A pair of permutations {P1_{m,k},P2_{m,k}} with P1[i] != P2[i]; i= 1..m1  
When the condition P1[i] != P2[i]; i= 1..m1 is satisfied the resulting latin squares are automatically orthogonal themselves, such that no double numbers occur in the resulting square 

these pairs are yet to be found, with the linear it was simple selecting an other parameter, this new situation needs some further study; 

The amount of orthogonal permutation is the set of permutations the satisfy the pandiagonality conditon determine the amount of the "most basic" pandiagonal magic squares, another matter took too much time to determine that number. The argumenttion further remains the same as in the table above,some (near future (I hope)) upload more details will be given. 
prime order panmagic multiplicative squares  

panmagic multiplicative squares in relation to panmagic addititive squares  
panmagic additive square  S = [m * LS(a)_{=P(m)} + LS(b)_{=P(m)} + 1];a < m/2; b != a;  
See derivation in previous table  
panmagic multiplicative square  P = [p1^{(LS(a)=P(m))} * p2^{(LS(b)=P(m))}]; p1 != p2 (relatively prime)  
given the fact that the panmagic additive squares of prime order (given by S) consists of all the numbers in the range [1 .. n^{2}] and the fact that all composite numbers have a unique decomposition in primes ensures that all the numbers in P are unique numbers albeit with many a gap. 

magic product 
p = _{i=0}∏^{m1}P_{i,j} =
_{i=0}∏^{m1} [p1^{LS(a)} * p2^{LS(b)}] =
[p1^(_{i=0}∑^{m1} LS(a)) * p2^(_{i=0}∑^{m1}LS(b))] = p1^{m(m1)/2} p2^{m(m1)/2} = (p1 p2) ^{m(m1)/2} 

Although the above is worked out in only 1 1agonal direction the derivation hold on any 1agonal and (broken) 2agonal the crucial step is the fact that _{i=0}∑^{m1} LS(a) on each line sums over all the digits [0 .. m1] on any line, henche this sums to m(m1)/2 Seen the construction of these squares the amount of square is the same of the amount of additive magic squares for every chosen pair of relative primes p1 and p2, below the magic products are calculated for the case p1 = 2 and p2 = 3 

(2 * 3) ^{m(m1)/2}  
5 7 11 13 
60.466.176 21.936.950.640.377.900 6.285.195.213.566 * 10^{30} 4.963.608.617.944.920 * 10^{45} 