import numpy as np
    import matplotlib.pyplot as plt
    import glob

    def get_data(fname):
        with open(fname, 'r') as infile:
            lines = infile.readlines()
        SOE = [i for (i, line) in enumerate(lines) if 'SOE' in line][0]
        EOE = [i for (i, line) in enumerate(lines) if 'EOE' in line][0]
        lines = lines[SOE+1:EOE]
        data = np.array([[float(x) for x in line.split(',')[2:8]]
                         for line in lines])
        position, velocity = data.T.reshape(2, 3, -1)
        distance, speed = [np.sqrt((thing**2).sum(axis=0))
                           for thing in (position, velocity)]
        JD = np.array([float(line.split(',')[0]) for line in lines])
        return lines, JD, distance, speed

    fnames = sorted(glob.glob('* horizons_results.txt'))

    titles = 'SS Barycenter', 'Sun'

    fig, (ax_bary, ax_Sun) = plt.subplots(2, 1)

    d_bary, d_Sun = [], []


    for fname in fnames:
        name, observer = fname.split()[:2]
        lines, JD, distance, speed = get_data(fname)
        if 'barycenter' in observer.lower():
            ax_bary.plot(JD-JD[0]+1, distance/1E+09, label=name)
            d_bary.append(distance)
        else:
            ax_Sun.plot(JD-JD[0]+1, distance/1E+09, label=name)
            d_Sun.append(distance)
        ax_bary.set_title('SS Barycenter')
        ax_bary.set_ylabel('distance (km×10⁹)')
        ax_bary.legend()
        ax_Sun.set_title('Sun')
        ax_Sun.set_ylabel('distance (km×10⁹)')
        ax_Sun.legend()
        ax_Sun.set_xlabel('day in July 2023 (UTC)')
    diff_bary = d_bary[1] - d_bary[0]
    diff_Sun = d_Sun[1] - d_Sun[0]
    i_bary = np.argmax(diff_bary > 0)
    i_Sun = np.argmax(diff_Sun > 0)
    ax_bary.plot(JD[i_bary:i_bary+1]-JD[0]+1, d_bary[0][i_bary:i_bary+1]/1E+09, 'ok')
    ax_Sun.plot(JD[i_Sun:i_Sun+1]-JD[0]+1, d_Sun[0][i_Sun:i_Sun+1]/1E+09, 'ok')
    xy = (JD[i_bary]-JD[0]+1, d_bary[0][i_bary]/1E+09)
    print('bary: ', xy[1])
    date = lines[i_bary].split(',')[1][:-8]
    ax_bary.annotate(date, xy)
    xy = (JD[i_Sun]-JD[0]+1, d_Sun[0][i_Sun]/1E+09)
    print('Sun: ', xy[1])
    date = lines[i_Sun].split(',')[1][:-8]
    ax_Sun.annotate(date, xy)
    plt.show()