#! /usr/bin/env python3
#
def problem1 ( ):

#*****************************************************************************80
#
## problem1(): Verify A*A is the identity matrix.
# 
  import numpy as np
 
  print ( '' )
  print ( 'problem1():' )
  print ( '  Define matrix A, verify that A*A = I' )
  print ( '' )

  n = 5
  A = problem1_matrix ( n )

  I = np.matmul ( A, A )

  print ( '  A*A = ', I )

  return

def problem1_matrix ( n ):

  import numpy as np

  A = np.zeros ( [ n, n ] )
  for i in range ( 0, n ):
    for j in range ( 0, n ):
      A[i,j] = np.sqrt ( 2 / ( n + 1 ) ) \
        * np.sin ( ( i + 1 ) * ( j + 1 ) * np.pi / ( n + 1 ) )

  return A

  return
  
def problem2 ( ):

#*****************************************************************************80
#
## problem2(): Draw a hexagonal figure divided into colored triangles.
#
  import matplotlib.pyplot as plt
  
  print ( '' )
  print ( 'problem2():' )
  print ( '  Draw a hexagonal figure divided into colored triangles.' )
  print ( '' )

  ax = -1.0; ay = 0.0;
  bx = 1.0; by = 0.0;
  cx = 2.0; cy = 1.0;
  dx = 1.0; dy = 2.0;
  ex = -1.0; ey = 2.0;
  fx = -2.0; fy = 1.0;
  zx = 0.0; zy = 1.0;

  plt.clf ( )
  plt.fill ( [ ax, bx, zx ], [ ay, by, zy ], 'r' )
  plt.fill ( [ bx, cx, zx ], [ by, cy, zy ], 'g' )
  plt.fill ( [ cx, dx, zx ], [ cy, dy, zy ], 'b' )
  plt.fill ( [ dx, ex, zx ], [ dy, ey, zy ], 'r' )
  plt.fill ( [ ex, fx, zx ], [ ey, fy, zy ], 'g' )
  plt.fill ( [ fx, ax, zx ], [ fy, ay, zy ], 'b' )

  plt.grid ( True )
  plt.axis ( 'equal' )
  filename = 'problem2.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.close ( )

  return

def problem3 ( ):

#*****************************************************************************80
#
## problem3(): Draw a triangular grid of blue dots.
#  
  import matplotlib.pyplot as plt
  
  print ( '' )
  print ( 'problem3():' )
  print ( '  Draw a triangular grid of blue dots.' )
  print ( '' )

  plt.clf ( )

  for i in range ( 0, 11 ):
    for j in range ( 0, i + 1 ):
      plt.plot ( i, j, 'b.', markersize = 20 )

  plt.grid ( True )
  plt.axis ( 'equal' )
  filename = 'problem3.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.close ( )

  return
  
def problem4 ( ):

#*****************************************************************************80
#
## problem4(): plot the sinc function, with a "cutoff" line.
#
  from math import sin
  import matplotlib.pyplot as plt

  print ( '' )
  print ( 'problem4():' )
  print ( '  plot the sinc function, with a "cutoff" line..' )
  print ( '' )

  n = 101
  xmin = -10.0
  xmax = +10.0
  xlist = []
  ylist = []
  for i in range ( 0, 101 ):
    x = ( ( n - i ) * xmin + i * xmax ) / ( n - 1 )
    y = sin ( x ) / x
    xlist.append ( x )
    ylist.append ( y )

  plt.clf ( )
  plt.plot ( xlist, ylist )
  plt.plot ( [xmin,xmax], [0.8,0.8], 'r-', linewidth = 2 )
  plt.grid ( True )
  filename = 'problem4.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.close ( )

  return

def problem5 ( ):

#*****************************************************************************80
#
## problem5(): complex arithmetic
#
  from math import atan2
  from math import cos
  from math import sin

  print ( '' )
  print ( 'problem5():' )
  print ( '  Complex arithmetic' )
  print ( '' )

  ax = -5
  ay = 4
  a = complex ( ax, ay )

  print ( '  a = ', a )

  bx = 2
  by = -3
  b = complex ( bx, by )
  print ( '  b = ', b )

  c = a * b
  cx = c.real
  cy = c.imag
  print ( '  c = ', c )
#
#  R:
#
  ar = abs ( a )
  br = abs ( b )
  cr = abs ( c )
  print ( '  ar = ', ar, '  br = ', br, '  cr = ', cr )
#
#  Verify cr = ar * br
#
  print ( '  cr = ', cr, '  ar * br = ', ar * br )
#
#  Theta
#
  at = atan2 ( ay, ax )
  bt = atan2 ( by, bx )
  ct = atan2 ( cy, cx )

  print ( '  at = ', at, '  bt = ', bt, '  ct = ', ct )
#
#  Verify ct = at + bt mod 2pi
#
  print ( '  ct = ', ct, '  at + bt = ', at + bt )
#
#  Polar
#
  ap = ar * complex ( cos ( at ), sin ( at ) )
  bp = br * complex ( cos ( bt ), sin ( bt ) )
  cp = cr * complex ( cos ( ct ), sin ( ct ) )
  print ( '  ap = ', ap, '  bp = ', bp, '  cp = ', cp )

  return
 
def problem6 ( ):

#*****************************************************************************80
#
## problem6(): minimum and maximum of a function.
#
  import numpy as np

  print ( '' )
  print ( 'problem6():' )
  print ( '  Find minimum and maximum of f(x) over an interval.' )
  print ( '' )

  xmin = 0.0
  xmax = 1.0
  n = 101
  x = np.linspace ( xmin, xmax, n )
  y = problem6_f ( x )

  miny = np.min ( y )
  maxy = np.max ( y )
  mini = np.argmin ( y )
  maxi = np.argmax ( y )
  minx = x[mini]
  maxx = x[maxi]

  print ( '  minimum occurs at x = ', x[mini], ' where y = ', miny )
  print ( '  maximum occurs at x = ', x[maxi], ' where y = ', maxy )

  return

def problem6_f ( x ):

  import numpy as np

  value = np.cos (  7.0 * x ) \
    + 5 * np.cos ( 11.2 * x ) \
    - 2 * np.cos ( 14.0 * x ) \
    + 5 * np.cos ( 31.5 * x ) \
    + 7 * np.cos ( 63.0 * x )

  return value
  
def problem7 ( ):

#*****************************************************************************80
#
## problem7(): Is a vector monotonically increasing?
#
  import numpy as np

  print ( '' )
  print ( 'problem7():' )
  print ( '  Write a one line statement that is:' )
  print ( '  * True if vector x is monotonically increasing. ' )
  print ( '  * True if vector x is strictly monotonically increasing.' )
  print ( '' )
  print ( '  np.diff(x) returns the successive difference of entries in x' )
  print ( '  What condition must this difference vector satisfy?' )
  print ( '  Hill Q6.1.10, page 226.' )
  print ( '' )
  print ( '  Test on the following:' )
  print ( '  1) x = [  1, 2, 3, 6, 8 ]' )
  print ( '  2) y = [ -1, 2, 3, 3, 7 ]' )
  print ( '  3) z = [  1, 3, 7, 2, 9 ]' )
 
  x = np.array ( [  1, 2, 3, 6, 8 ] )
  y = np.array ( [ -1, 2, 3, 3, 7 ] )
  z = np.array ( [  1, 3, 7, 2, 9 ] )
  xd = np.diff ( x )
  yd = np.diff ( y )
  zd = np.diff ( z )

  print ( '  x monotonic          is ', np.all ( 0 <= xd ) )
  print ( '  x strictly monotonic is ', np.all ( 0 <  xd ) )
  print ( '  y monotonic          is ', np.all ( 0 <= yd ) )
  print ( '  y strictly monotonic is ', np.all ( 0 <  yd ) )
  print ( '  z monotonic          is ', np.all ( 0 <= zd ) )
  print ( '  z strictly monotonic is ', np.all ( 0 <  zd ) )

  return
  
def problem8 ( ):

#*****************************************************************************80
#
## problem8(): Area of a polygon.
#
  print ( '' )
  print ( 'problem8():' )
  print ( '  Area of a polygon (Hill P6.1.2.' )
  print ( '  Use for() loops and a list.' )
  print ( '' )

  n = 8
  x = [ 0, 3, 3, 2, 2, 1, 1, 0 ]
  y = [ 0, 0, 3, 3, 1, 1, 2, 2 ]

  x.append ( x[0] )
  y.append ( y[0] )

  s1 = 0
  for i in range ( 0, n ):
    s1 = s1 + x[i] * y[i+1]

  s2 = 0
  for i in range ( 0, n ):
    s2 = s2 + x[i+1] * y[i]

  area = ( s1 - s2 ) / 2
  print ( '  Area = ', area )

  return

def problem9 ( ):

#*****************************************************************************80
#
## problem9(): Area of a polygon.
#
  import numpy as np

  print ( '' )
  print ( 'problem9():' )
  print ( '  Use numpy arrays and indexing.' )
  print ( '' )

  n = 8
  x = np.array ( [ 0, 3, 3, 2, 2, 1, 1, 0 ] )
  y = np.array ( [ 0, 0, 3, 3, 1, 1, 2, 2 ] )

  np.append ( x, x[0] )
  np.append ( y, y[0] )

  area = 0.5 * ( np.sum ( x[0:n-1] * y[1:n] ) - np.sum ( x[1:n] * y[0:n-1] ) )
  print ( '  Area = ', area )

  return

if ( __name__ == "__main__" ):
  problem1 ( )
  problem2 ( )
  problem3 ( )
  problem4 ( )
  problem5 ( )
  problem6 ( )
  problem7 ( )
  problem8 ( )
  problem9 ( )