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.