#! /usr/bin/env python3
#
def arenstorf_dydt ( t, xy ):

#*****************************************************************************80
#
## arenstorf_dydt() evaluates the right hand side of the Arenstorf ODE.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    07 November 2024
#
#  Author:
#
#    John Burkardt
#
  import numpy as np

  global mu1, mu2

  x = xy[0]
  y = xy[1]
  xp = xy[2]
  yp = xy[3]

  d1 = np.sqrt ( ( ( x + mu1 )**2 + y**2 )**3 )
  d2 = np.sqrt ( ( ( x - mu2 )**2 + y**2 )**3 )

  dxdt = xp
  dydt = yp
  dxpdt = x + 2.0 * yp - mu2 * ( x + mu1 ) / d1 - mu1 * ( x - mu2 ) / d2
  dypdt = y - 2.0 * xp - mu2 * y / d1 - mu1 * y / d2

  dxydt = np.array ( [ dxdt, dydt, dxpdt, dypdt ] )

  return dxydt

def arenstorf_solve_ivp ( ):

#*****************************************************************************80
#
## arenstorf_solve_ivp() uses solve_ivp() to solve the Arenstorf ODE.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    07 November 2024
#
#  Author:
#
#    John Burkardt
#
  from scipy.integrate import solve_ivp
  import matplotlib.pyplot as plt
  import numpy as np

  global mu1, mu2

  mu1 = 0.012277471
  mu2 = 1.0 - mu1
  tmin = 0.0
  tmax = 17.0652165601579625588917206249
  n = 71

  tspan = np.array ( [ tmin, tmax ] )
  t = np.linspace ( tmin, tmax, n )
  y0 = np.array ( [ 0.994, 0.0, 0.0, -2.00158510637908252240537862224 ] )

  sol = solve_ivp ( arenstorf_dydt, tspan, y0, t_eval = t )
#
#  Set the moon's location.
#
  theta = np.linspace ( 0.0, 2.0 * np.pi, n )
  moon_x = np.cos ( theta )
  moon_y = np.sin ( theta )
#
#  Plot the orbit as a sequence of frames.
#
  for i in range ( 0, n ):
#
#  Plot earth.
#
    plt.plot ( 0.0, 0.0, 'bo', markersize = 25 )
#
#  Plot moon.
#
    plt.plot ( moon_x[0:i+1], moon_y[0:i+1], 'go', markersize = 15 )
#
#  Plot satellite.
#
    plt.plot ( sol.y[0,0:i+1], sol.y[1,0:i+1], 'ro', markersize = 10 )
#
#  Finish the plot.
#
    plt.grid ( True )
    plt.axis ( 'Equal' )
    plt.xlabel ( '<---  x  --->' )
    plt.ylabel ( '<---  y  --->' )
    plt.title ( 'Arenstorf orbit' )
    filename = 'arenstorf_solve_ivp_phase.png'
    plt.savefig ( filename )
    print ( '  Graphics saved as "%s"' % ( filename ) )
    plt.show ( )
    plt.close ( )

  return

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the MIT license. 
#
#  Modified:
#
#    21 August 2019
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return

if ( __name__ == '__main__' ):
  timestamp ( )
  arenstorf_solve_ivp ( )
  timestamp ( )