#! /usr/bin/env python3
#
def magic_matrix ( n ):

#*****************************************************************************80
#
## magic_matrix() returns a magic matrix of odd order n.
#
#  Discussion:
#
#    Every row and column of the matrix has the same sum.
#
#    The algorithm proceeds as follows:
#
#    a) Start in the middle of the top row, and let k = 1
#    b) Insert k into the grid.
#    c) If k = n^2, you are done.
#    d) Increment k.
#    e) Move diagonally up and to the right, but wrap to the first
#       column, or to the last row, if you leave the grid.
#    f) If that cell is already occupied, drop down one space from
#       your current position.
#    g) Return to step (b).  
# 
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    05 December 2022
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    integer n: the number of rows and columns in the matrix,
#    which must be odd.
#
#  Output:
#
#    integer A(n,n): the magic matrix.
#
  import numpy as np

  if ( ( n % 2 ) != 1 ):
    print ( '' )
    print ( 'magic_matrix(): Fatal error!' )
    print ( '  Input value n must be an odd integer.' )
    raise Exception ( 'magic_matrix(): Fatal error!' )

  A = np.zeros ( [ n, n ], dtype = int )

  k = 1
  i = 1
  j = ( n + 1 ) // 2
  A[i-1,j-1] = k

  while ( k < n**2 ):

    k = k + 1

    im1 = ( ( i - 2 ) % n ) + 1
    jp1 = ( j % n ) + 1
    if ( A[im1-1,jp1-1] != 0 ):
      im1 = i + 1
      jp1 = j

    A[im1-1,jp1-1] = k
    i = im1
    j = jp1

  return A

def magic_matrix_test ( ):

#*****************************************************************************80
#
## magic_matrix_test() tests magic_matrix().
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    05 December 2022
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'magic_matrix_test():' )
  print ( '  Python version: ' + platform.python_version ( ) )
  print ( '  Test magic_matrix()' )

  for n in [ 3, 5 ]:
    A = magic_matrix ( n )
    print ( '' )
    print ( '  Magic matrix for n =', n )
    print ( A )
    print ( '  Row and column sums:' )
    rsum = np.sum ( A, axis = 1 )
    print ( rsum )
    csum = np.sum ( A, axis = 0 )
    print ( csum )
#
#  Terminate.
#
  print ( '' )
  print ( 'magic_matrix_test():' )
  print ( '  Normal end of execution.' )

  return

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    21 August 2019
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return

if ( __name__ == '__main__' ):
  timestamp ( )
  magic_matrix_test ( )
  timestamp ( )