#! /usr/bin/env python3
#
def tpy_direction_field ( show_solution ):

  import matplotlib.pyplot as plt
  import numpy as np
#
#  Create a -10 x 10 TY grid with unit spacing.
#
  tmin = -10.0
  tmax = +10.0
  ymin = -10.0
  ymax = +10.0

  tvec = np.linspace ( tmin, tmax, 21 )
  dt = tvec[1] - tvec[0]
  yvec = np.linspace ( ymin, ymax, 21 )
  T, Y = np.meshgrid ( tvec, yvec )
#
#  Define DT, DY arrays and plot the direction field.
#
  DT = dt * np.ones ( T.shape )
  DY = dt * ( 0.2 * T + 0.125 * Y )
#
#  Make this true to plot normalized direction vectors.
#
  if ( False ):
    S = np.sqrt ( DT**2 + DY**2 )
    DT = DT / S
    DY = DY / S
#
#  Prepare for plotting.
#
  plt.clf ( )
#
#  Plot the flow field or direction field.
#
  plt.quiver ( T, Y, DT, DY, color = 'red' )
#
#  Limit plotting to the [tmin,tmax] x [ymin,ymax] box.
#
  plt.xlim ( [ tmin, tmax ] )
  plt.ylim ( [ ymin, ymax ] )
#
#  Draw axis lines.
#
  plt.plot ( [ 0.0, 0.0 ], [ ymin, ymax ], 'k-' )
  plt.plot ( [ tmin, tmax ], [ 0.0, 0.0 ], 'k-' )
#
#  Use the Euler method to estimate the solution 
#  for several initial values y0.
#
#  The equation is dy/dt = 1/5 t + 1/8 y 
#
  if ( show_solution ):
    n_euler = 21
    t = np.linspace ( tmin, tmax, n_euler )
    dt = t[1] - t[0]
    y = np.zeros ( n_euler )

    for y0 in range ( 0, 11 ):
      for i in range ( 0, n_euler ):
        if ( i == 0 ):
          y[i] = y0
        else:
          y[i] = y[i-1] + dt * ( 0.2 * t[i-1] + 0.125 * y[i-1] )
      plt.plot ( t, y, linewidth = 3 )

  plt.grid ( True )
  if ( show_solution ):
    filename = 'tpy_solution_field.png'
  else:
    filename = 'tpy_direction_field.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( )
  plt.close ( )

  return

if ( __name__ == "__main__" ):
  tpy_direction_field ( False )
  tpy_direction_field ( True )