simple planetary simulation

planet.py

import numpy as np

def clip_norm(x):
    return x / abs(np.max(x)) * 0.001

# random strat
# def p1_strat(p, v):
#     fr = np.random.uniform(-1.0, 1.0, size=(2,))
#     return fr
def p2_strat(p, v):
    fr = np.random.uniform(-1.0, 1.0, size=(2,))
    return fr

# projective strat that goes toward goal
def p1_strat(p, v):
    delta1 = np.array([-1.0, 0.0]) - p
    f1 = 0.001 / (0.1 + np.sum(delta1**2)) * delta1
    delta2 = np.array([1.0, 0.0]) - p
    f2 = 0.001 / (0.1 + np.sum(delta2**2)) * delta2
    goal = np.array([-1.0, 0.0]) - p
    to_goal = goal - (f1 + f2 + v)
    to_goal_max = abs(np.max(to_goal))
    return to_goal / to_goal_max * 0.001


# 2 step look ahead of looking at p1's strat
# def p2_strat(p, v):
#     delta1 = np.array([-1.0, 0.0]) - p
#     f1 = 0.001 / (0.1 + np.sum(delta1**2)) * delta1
#     delta2 = np.array([1.0, 0.0]) - p
#     f2 = 0.001 / (0.1 + np.sum(delta2**2)) * delta2
# 
#     v_p1 = p1_strat(p, v)
# 
#     goal = np.array([1.0, 0.0]) - p
#     to_goal = goal - (f1 + f2 + v + v_p1)
#     to_goal_max = abs(np.max(to_goal))
#     ret_greedy = to_goal / to_goal_max * 0.001
# 
#     # subversion
#     if np.dot((f1 + f2 + v + v_p1), delta1) > 0 and np.sum(delta1 ** 2) < 0.4:
#         # v_max = abs(np.max(v)) + 1e-5
#         # return v / v_max * 0.001
#         cur_dir = (f1 + f2 + v + v_p1)
#         per1 = np.array([cur_dir[1], -cur_dir[0]])
#         per2 = np.array([-cur_dir[1], cur_dir[0]])
#         if np.dot(per1, v_p1) > np.dot(per2, v_p1):
#             return clip_norm(per2)
#         else:
#             return clip_norm(per1)
#     # attraction
#     else:
#         return ret_greedy

def step(pts, vs):
    new_ps = []
    new_vs = []
    for i in range(len(pts)):
        p = pts[i]
        v = vs[i]
        new_pp = p + v
        # mess with the velocity to tend to ward the center
        
        # force to 2 attractors:
        delta1 = np.array([-1.0, 0.0]) - p
        f1 = 0.001 / (0.1 + np.sum(delta1**2)) * delta1
        delta2 = np.array([1.0, 0.0]) - p
        f2 = 0.001 / (0.1 + np.sum(delta2**2)) * delta2

        # a random force representing the action of the players

        move1 = np.clip(p1_strat(p, v), -0.001, 0.001)
        move2 = np.clip(p2_strat(p, v), -0.001, 0.001)

        new_v = v + f1 + f2 + move1 + move2
        new_v = 0.99 * new_v

        if np.sum(delta1 ** 2) < 0.005 or np.sum(delta2 ** 2) < 0.005:
            new_v = 0.0 * new_v
        new_vs.append(new_v)
        new_ps.append(new_pp)
    return new_ps, new_vs

ps = [np.random.uniform(-2.0, 2.0, size=(2,)) for _ in range(100)]
vs = [0.1 * np.random.uniform(-1.0, 1.0, size=(2,)) for _ in range(100)]

# ps = [np.array([-2 + 0.1 * i, 1.9]) for i in range(40)]
# vs = [np.array([0.0, -0.01]) for _ in range(40)]

import matplotlib.pyplot as plt

for i in range(1200):
    print (i)
    axes = plt.gca()
    axes.set_xlim([-2 , 2])
    axes.set_ylim([-2 , 2])
    plt.scatter([x[0] for x in ps], [x[1] for x in ps])
    plt.scatter([-1, 1], [0, 0], c='red')
    plt.savefig('figs/'+str(i)+".png")
    plt.clf()
    ps, vs = step(ps, vs)