#! /usr/bin/env python
#
def moment_central ( n, x, y, p, q ):

#*****************************************************************************80
#
## MOMENT_CENTRAL computes central moments of a polygon.
#
#  Discussion:
#
#    The central moment Mu(P,Q) is defined by
#
#      Mu(P,Q) = Integral ( polygon ) (x-Alpha(1,0))^p (y-Alpha(0,1))^q dx dy
#              / Area ( polygon )
#
#    where
#
#      Alpha(1,0) = Integral ( polygon ) x dx dy / Area ( polygon )
#      Alpha(0,1) = Integral ( polygon ) y dx dy / Area ( polygon )
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 June 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Carsten Steger,
#    On the calculation of arbitrary moments of polygons,
#    Technical Report FGBV-96-05,
#    Forschungsgruppe Bildverstehen, Informatik IX,
#    Technische Universitaet Muenchen, October 1996.
#
#  Parameters:
#
#    Input, integer N, the number of vertices of the polygon.
#
#    Input, real X(N), Y(N), the vertex coordinates.
#
#    Input, integer P, Q, the indices of the moment.
#
#    Output, real MU_PQ, the uncentral moment Mu(P,Q).
#
  from moment_normalized import moment_normalized
  from r8_choose import r8_choose
  from r8_mop import r8_mop

  alpha_10 = moment_normalized ( n, x, y, 1, 0 )
  alpha_01 = moment_normalized ( n, x, y, 0, 1 )

  mu_pq = 0.0

  for i in range ( 0, p + 1 ):
    for j in range ( 0, q + 1 ):

      alpha_ij = moment_normalized ( n, x, y, i, j )

      mu_pq = mu_pq + r8_mop ( p + q - i - j ) \
        * r8_choose ( p, i ) * r8_choose ( q, j ) \
        * alpha_10 ** ( p - i ) * alpha_01 ** ( q - j ) * alpha_ij

  return mu_pq

def moment_central_test ( ):

#*****************************************************************************80
#
## MOMENT_CENTRAL_TEST tests MOMENT_CENTRAL.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 June 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  n = 4

  mu_exact = np.array ( [ \
    1.0, 0.0, 0.0, 5.666666666666667, 2.0, 2.666666666666667 ] )

  x = np.array ( [ 2.0, 10.0, 8.0, 0.0 ] )
  y = np.array ( [ 0.0,  4.0, 8.0, 4.0 ] )

  print ( '' )
  print ( 'MOMENT_CENTRAL_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  MOMENT_CENTRAL computes central moments of a polygon.' )
  print ( '  Here, we test the code using a rectange with known moments.' )
  print ( '' )
  print ( '   P   Q             Mu(P,Q)' )
  print ( '            Computed         Exact' )
  print ( '' )
  k = 0
  for s in range ( 0, 3 ):
    for p in range ( s, -1, -1 ):
      q = s - p
      mu_pq = moment_central ( n, x, y, p, q )
      print ( '  %2d  %2d  %14.6g  %14.6g' % ( p, q, mu_pq, mu_exact[k] ) )
      k = k + 1
#
#  Terminate.
#
  print ( '' )
  print ( 'MOMENT_CENTRAL_TEST' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  moment_central_test ( )
  timestamp ( )