Ed Pegg:
One can use a modified De LaLoubere (Siam) algorithm to reach a similar pandiagonal magic square. Where the matrix size is 7, and place first number in a normal sequence at cell 2,3. Then apply the following rules to fill the matrix using ordinary step 1 over to the right with ordinary step 4 down. If landing cell is not empty, then make a break step 5 over to right with break step 4 up place next number in a normal sequence.
The following code is old and never really modified it to flow better in WL, it is a translation from a Pascal version I wrote >25 years ago. So does not represent efficient WL code.
magicSquare[n_, ordx_, ordy_, breakx_, breaky_, i_, j_] :=
Module[{magicMatrix = Table[0, {n}, {n}], x = i, y = j, counter = 1,
nextox, nextoy, nextbx, nextby}, magicMatrix[[x, y]] = 1;
While[counter < n^2,(*while body*)counter++;
nextox = x + ordy;
If[nextox < 1, nextox = n + nextox,];
If[nextox > n, nextox = Mod[nextox, n],];
nextoy = y + ordx;
If[nextoy < 1, nextoy = n + nextoy,];
If[nextoy > n, nextoy = Mod[nextoy, n],];
If[magicMatrix[[nextox, nextoy]] == 0,(*if body*)
x = nextox; y = nextoy;
magicMatrix[[x, y]] = counter,
(*else*)nextbx = x + breaky;
If[nextbx < 1, nextbx = n + nextbx,];
If[nextbx > n, nextbx = Mod[nextbx, n],];
nextby = y + breakx;
If[nextby < 1, nextby = n + nextby,];
If[nextby > n, nextby = Mod[nextby, n],];
If[magicMatrix[[nextbx, nextby]] == 0, x = nextbx;
y = nextby;
magicMatrix[[x, y]] = counter],];(*end if*)];(*end while*)
Return[MatrixForm[magicMatrix]];];
So applying the function
z = magicSquare[7, 1, 4, 5, -4, 2, 3]
We get
(6 46 37 35 26 17 8
19 10 1 48 39 30 28
32 23 21 12 3 43 41
45 36 34 25 16 14 5
9 7 47 38 29 27 18
22 20 11 2 49 40 31
42 33 24 15 13 4 44 )
Which seems equivalent to the matrix in your example. If you Mod[ z[[[1]],7] +1, then you get a Latin Square.