kepler_interpolate_sample

A sample script for Kepler class.

''' A sample script for Kepler class.
'''

import os
import sys
import logging
logging.basicConfig()
import datetime
import math

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

import irfpy.util.keplernumeric
import irfpy.util.ballistic

mjup = 1.8986e27
mio = 8.93e22

rio0 = 421700e3
vio0 = 17e3

delta_t = 1080
total_t = 60 * 3600
n_proc = int(total_t / delta_t)

def main():
    # Instance Kepler.
    kepler = irfpy.util.keplernumeric.Kepler(mjup, mio, np.array([rio0, 0, 0]), np.array([0, vio0, 0]), delta_t)  # Every 0.3 hours.

    kepler.proceed(n_proc)   # In total 60 hours.

    tx, xyz = kepler.get_xlist()

    plt.plot(xyz[:, 0], xyz[:, 1], 'x', label='Leapflog')

    # Interpolation
    tlist = np.linspace(0, 86400, 86400)
    xyzinter = kepler.get_pos(tlist)

    plt.plot(xyzinter[0, :], xyzinter[1, :], '-', label='LF/INTP')

    # Backtarcing
    kepler_bs = irfpy.util.keplernumeric.Kepler(mjup, mio, np.array([rio0, 0, 0]), np.array([0, vio0, 0]), -delta_t)
    kepler_bs.proceed(n_proc // 4)    # In total 15 hours
    tlist = np.linspace(-25000, 0, 5000)
    xyzinter = kepler_bs.get_pos(tlist)

    plt.plot(xyzinter[0, :], xyzinter[1, :], '-', label='LF/INTP (BS)')
    
    # Theoretical (ballistic)
    ballist = irfpy.util.ballistic.Trajectory(rio0, 0, 0, vio0 / rio0, mjup, m=mio)
    ballist.print_values()
    ballang = np.linspace(0, np.pi, 100)
    ballr = [ballist.r(t) for t in ballang]

    plt.plot(ballr * np.cos(ballang), ballr * np.sin(ballang), 'k-', lw=2, label='Ballist')

    plt.legend()
    plt.gca().set_aspect('equal')

    plt.axes([0.1, 0.1, 0.4, 0.2])

    #### Time series plot for numeric solver
    tlist = np.linspace(0, 86400, 86400)
    vxyzinter = kepler.get_vel(tlist)
    plt.plot(tlist, vxyzinter[0, :], 'b', label='Vx')
    plt.plot(tlist, vxyzinter[1, :], 'g', label='Vy')

    tlist = np.linspace(-25000, 0, 5000)
    vxyzinter = kepler_bs.get_vel(tlist)
    plt.plot(tlist, vxyzinter[0, :], 'b')
    plt.plot(tlist, vxyzinter[1, :], 'g')

    plt.legend()
    plt.savefig('kepler_interpolate_sample.png')


if __name__ == "__main__":
    main()