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

#*****************************************************************************80
#
## pendulum_solve_ivp() solves the pendulum ODE using solve_ivp().
#
  from scipy.integrate import solve_ivp
  from pendulum_conserved import pendulum_conserved
  from pendulum_dydt import pendulum_dydt
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'pendulum_solve_ivp:' )
  print ( '  Solve the pendulum ODE using solve_ivp()' )

  t0 = 0.0
  tstop = 20.0
  tspan = np.array ( [ t0, tstop ] )
  y0 = np.array ( [ np.pi / 3.0, 0.0 ] )

  sol = solve_ivp ( pendulum_dydt, tspan, y0 )

  y = np.transpose ( sol.y )
#
#  Time plot.
#
  plt.clf ( )
  plt.plot ( sol.t, y, linewidth = 3 )
  plt.grid ( True );
  plt.xlabel ( '<-- Time -->' )
  plt.ylabel ( '<-- Theta, dThetadt -->' )
  plt.legend ( [ 'Theta', 'dThetadt' ] )
  plt.title ( 'pendulum, solve_ivp, time plot' )
  filename = 'pendulum_solve_ivp_time.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( )
  plt.close ( )
#
#  Conservation plot.
#
  h = pendulum_conserved ( sol.t, y )

  plt.clf ( )
  plt.plot ( sol.t, h, 'c-', lineWidth = 3 )
  plt.plot ( [ t0, tstop ], [ 0.0, 0.0 ], 'k-', linewidth = 3 )
  plt.grid ( 'on' )
  plt.xlabel ( '<-- Time -->' )
  plt.ylabel ( '<-- Energy(t) -->' )
  plt.title ( 'pendulum, solve_ivp, conservation plot' )
  filename = 'pendulum_solve_ivp_conservation.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( )
  plt.close ( )

  return

if ( __name__ == "__main__" ):
  pendulum_solve_ivp ( )