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

#*****************************************************************************80
#
## MOMENT computes an unnormalized moment of a polygon.
#
#  Discussion:
#
#    Nu(P,Q) = Integral ( x, y in polygon ) x^p y^q dx dy
#
#  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 NU_PQ, the unnormalized moment Nu(P,Q).
#
  from r8_choose import r8_choose

  nu_pq = 0.0

  xj = x[n-1]
  yj = y[n-1]

  for i in range ( 0, n ):

    xi = x[i]
    yi = y[i]

    s_pq = 0.0

    for k in range ( 0, p + 1 ):
      for l in range ( 0, q + 1 ):
        s_pq = s_pq \
          + r8_choose ( k + l, l ) * r8_choose ( p + q - k - l, q - l ) \
          * xi ** k * xj ** ( p - k ) \
          * yi ** l * yj ** ( q - l )

    nu_pq = nu_pq + ( xj * yi - xi * yj ) * s_pq

    xj = xi
    yj = yi

  nu_pq = nu_pq / float ( p + q + 2 ) / float ( p + q + 1 ) \
    / r8_choose ( p + q, p )

  return nu_pq

def moment_test ( ):

#*****************************************************************************80
#
## MOMENT_TEST tests MOMENT.
#
#  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

  nu_exact = np.array ( [ \
    40.0, 200.0, 160.0, 1226.66666666666667, 880.0, 746.66666666666666 ] )

  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_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  MOMENT computes moments of a polygon.' )
  print ( '  Here, we test the code using a rectange with known moments.' )
  print ( '' )
  print ( '   P   Q             Nu(P,Q)' )
  print ( '            Computed         Exact' )
  print ( '' )
  k = 0
  for s in range ( 0, 3 ):
    for p in range ( s, -1, -1 ):
      q = s - p
      nu_pq = moment ( n, x, y, p, q )
      print ( '  %2d  %2d  %14.6g  %14.6g' % ( p, q, nu_pq, nu_exact[k] ) )
      k = k + 1
#
#  Terminate.
#
  print ( '' )
  print ( 'MOMENT_TEST' )
  print ( '  Normal end of execution.' )
  return

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