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

#*****************************************************************************80
#
## sphere_solve_ivp() solves the sphere ODE using solve_ivp().
#
  from scipy.integrate import solve_ivp
  from sphere_conserved import sphere_conserved
  from sphere_dydt import sphere_dydt
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'sphere_solve_ivp():' )
  print ( '  Solve the sphere ODE using solve_ivp()' )

  t0 = 0.0
  tstop = 45.0
  tspan = np.array ( [ t0, tstop ] )
  y0 = np.array ( [ np.cos ( 0.9 ), 0.0, np.sin ( 0.9 ) ] )

  t = np.linspace ( t0, tstop, 501 )

  sol = solve_ivp ( sphere_dydt, tspan, y0, t_eval = t )

  y = np.transpose ( sol.y )
#
#  Time plot.
#
  plt.clf ( )
  plt.plot ( sol.t, y, linewidth = 3 )
  plt.grid ( True );
  plt.xlabel ( '<-- Time -->' )
  plt.ylabel ( '<-- Position -->' )
  plt.legend ( [ 'P', 'Q', 'R' ] )
  plt.title ( 'sphere model, time plot' )
  filename = 'sphere_solve_ivp_time.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.show ( )
  plt.close ( )
#
#  Conservation plot.
#
  P = sphere_conserved ( sol.t, y )

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

  return

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