import matplotlib.pyplot as plt
import numpy as np

if __name__ == "__main__":
    fig = plt.figure(figsize=(15.2, 8))
    ax = fig.add_subplot(111)

    R = 10  # planet radius
    pp = [0, 0]  # planet position
    ps = [40, 0]  # satellite position
    k = 2750  # force scale

    dist = lambda p1, p2: np.sqrt( (p1[0]-p2[0])**2 + (p1[1]-p2[1])**2 )  # Euclidean distance
    irsq = lambda p: 1 / dist(p, ps)**2   # inverse square distance
    gmag = lambda p: k*irsq(p)           # satellite's gravity force magnitude
    gdir = lambda p: np.divide(np.subtract(ps, p), dist(ps, p))  # satellite's gravity force direction
    gvec = lambda p: gmag(p)*gdir(p)     # satellite's gravity force vector
    
    gp = gvec(pp)  # satellite's gravity vector at center of planet
    tvec = lambda p: gvec(p)-gp  # satellite's tide force vector

    theta = np.linspace(-np.pi, np.pi, 50)  # circle internal angle
    for ti in theta:
        pi = [-R * np.cos(ti), R * np.sin(ti)]  # evaluation point on the perimeter (planet's surface)
        ti = tvec(pi)  # tide at evaluation point
        gi = gvec(pi)  # gravity at evaluation point
         
        ax.arrow(pp[0], pp[1], gp[0], gp[1], head_width=0.5/2, head_length=0.7/2, fc='b', color='b')
        ax.arrow(pi[0], pi[1], gi[0], gi[1], head_width=0.5/2, head_length=0.7/2, fc='b', color='b')
        ax.arrow(pi[0], pi[1], ti[0], ti[1], head_width=0.5/1, head_length=0.7/1, fc='r', color='r')

    ax.plot(pp[0], pp[1], 'ok')
    #ax.annotate(r"+", xytext=(0,0), size=30, xy=(0,0), ha="center", va="center")
    ax.annotate(r"O", xytext=(-1,0), size=30, xy=(0,0), ha="center", va="center")
    ax.annotate(r"S", xy=(25,0), xytext=(15,0),
                size=30, va="center",
                arrowprops=dict(arrowstyle="fancy", fc='k'))
    ax.set_xlim(-13, 25)
    ax.set_ylim(-11, 11)
    ax.plot(R*np.cos(theta), R*np.sin(theta), '-k', linewidth=3.)
    ax.set_aspect('equal', 'box')
    ax.axis('off')
    plt.savefig('field tidal.svg', bbox_inches='tight', pad_inches=.15, transparent='true')