diff --git a/solutions/travis.yml b/.travis.yml similarity index 100% rename from solutions/travis.yml rename to .travis.yml diff --git a/Pipfile b/Pipfile new file mode 100644 index 0000000..aede842 --- /dev/null +++ b/Pipfile @@ -0,0 +1,13 @@ +[[source]] +name = "pypi" +url = "https://pypi.org/simple" +verify_ssl = true + +[dev-packages] + +[packages] +pytest = "*" +numpy = "*" + +[requires] +python_version = "3.7" diff --git a/Pipfile.lock b/Pipfile.lock new file mode 100644 index 0000000..9956e7f --- /dev/null +++ b/Pipfile.lock @@ -0,0 +1,101 @@ +{ + "_meta": { + "hash": { + "sha256": "8c364fbe4be6bc611b2c1423f61b9f0a5c5a023d80ba3cc114fe49a5b29a4b8e" + }, + "pipfile-spec": 6, + "requires": { + "python_version": "3.7" + }, + "sources": [ + { + "name": "pypi", + "url": "https://pypi.org/simple", + "verify_ssl": true + } + ] + }, + "default": { + "atomicwrites": { + "hashes": [ + "sha256:03472c30eb2c5d1ba9227e4c2ca66ab8287fbfbbda3888aa93dc2e28fc6811b4", + "sha256:75a9445bac02d8d058d5e1fe689654ba5a6556a1dfd8ce6ec55a0ed79866cfa6" + ], + "version": "==1.3.0" + }, + "attrs": { + "hashes": [ + "sha256:69c0dbf2ed392de1cb5ec704444b08a5ef81680a61cb899dc08127123af36a79", + "sha256:f0b870f674851ecbfbbbd364d6b5cbdff9dcedbc7f3f5e18a6891057f21fe399" + ], + "version": "==19.1.0" + }, + "more-itertools": { + "hashes": [ + "sha256:2112d2ca570bb7c3e53ea1a35cd5df42bb0fd10c45f0fb97178679c3c03d64c7", + "sha256:c3e4748ba1aad8dba30a4886b0b1a2004f9a863837b8654e7059eebf727afa5a" + ], + "markers": "python_version > '2.7'", + "version": "==7.0.0" + }, + "numpy": { + "hashes": [ + "sha256:0e2eed77804b2a6a88741f8fcac02c5499bba3953ec9c71e8b217fad4912c56c", + "sha256:1c666f04553ef70fda54adf097dbae7080645435fc273e2397f26bbf1d127bbb", + "sha256:1f46532afa7b2903bfb1b79becca2954c0a04389d19e03dc73f06b039048ac40", + "sha256:315fa1b1dfc16ae0f03f8fd1c55f23fd15368710f641d570236f3d78af55e340", + "sha256:3d5fcea4f5ed40c3280791d54da3ad2ecf896f4c87c877b113576b8280c59441", + "sha256:48241759b99d60aba63b0e590332c600fc4b46ad597c9b0a53f350b871ef0634", + "sha256:4b4f2924b36d857cf302aec369caac61e43500c17eeef0d7baacad1084c0ee84", + "sha256:54fe3b7ed9e7eb928bbc4318f954d133851865f062fa4bbb02ef8940bc67b5d2", + "sha256:5a8f021c70e6206c317974c93eaaf9bc2b56295b6b1cacccf88846e44a1f33fc", + "sha256:754a6be26d938e6ca91942804eb209307b73f806a1721176278a6038869a1686", + "sha256:771147e654e8b95eea1293174a94f34e2e77d5729ad44aefb62fbf8a79747a15", + "sha256:78a6f89da87eeb48014ec652a65c4ffde370c036d780a995edaeb121d3625621", + "sha256:7fde5c2a3a682a9e101e61d97696687ebdba47637611378b4127fe7e47fdf2bf", + "sha256:80d99399c97f646e873dd8ce87c38cfdbb668956bbc39bc1e6cac4b515bba2a0", + "sha256:88a72c1e45a0ae24d1f249a529d9f71fe82e6fa6a3fd61414b829396ec585900", + "sha256:a4f4460877a16ac73302a9c077ca545498d9fe64e6a81398d8e1a67e4695e3df", + "sha256:a61255a765b3ac73ee4b110b28fccfbf758c985677f526c2b4b39c48cc4b509d", + "sha256:ab4896a8c910b9a04c0142871d8800c76c8a2e5ff44763513e1dd9d9631ce897", + "sha256:abbd6b1c2ef6199f4b7ca9f818eb6b31f17b73a6110aadc4e4298c3f00fab24e", + "sha256:b16d88da290334e33ea992c56492326ea3b06233a00a1855414360b77ca72f26", + "sha256:b78a1defedb0e8f6ae1eb55fa6ac74ab42acc4569c3a2eacc2a407ee5d42ebcb", + "sha256:cfef82c43b8b29ca436560d51b2251d5117818a8d1fb74a8384a83c096745dad", + "sha256:d160e57731fcdec2beda807ebcabf39823c47e9409485b5a3a1db3a8c6ce763e" + ], + "index": "pypi", + "version": "==1.16.3" + }, + "pluggy": { + "hashes": [ + "sha256:19ecf9ce9db2fce065a7a0586e07cfb4ac8614fe96edf628a264b1c70116cf8f", + "sha256:84d306a647cc805219916e62aab89caa97a33a1dd8c342e87a37f91073cd4746" + ], + "version": "==0.9.0" + }, + "py": { + "hashes": [ + "sha256:64f65755aee5b381cea27766a3a147c3f15b9b6b9ac88676de66ba2ae36793fa", + "sha256:dc639b046a6e2cff5bbe40194ad65936d6ba360b52b3c3fe1d08a82dd50b5e53" + ], + "version": "==1.8.0" + }, + "pytest": { + "hashes": [ + "sha256:3773f4c235918987d51daf1db66d51c99fac654c81d6f2f709a046ab446d5e5d", + "sha256:b7802283b70ca24d7119b32915efa7c409982f59913c1a6c0640aacf118b95f5" + ], + "index": "pypi", + "version": "==4.4.1" + }, + "six": { + "hashes": [ + "sha256:3350809f0555b11f552448330d0b52d5f24c91a322ea4a15ef22629740f3761c", + "sha256:d16a0141ec1a18405cd4ce8b4613101da75da0e9a7aec5bdd4fa804d0e0eba73" + ], + "version": "==1.12.0" + } + }, + "develop": {} +} diff --git a/build/lib/pycontest/__init__.py b/build/lib/pycontest/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/build/lib/pycontest/elastic_collisions.py b/build/lib/pycontest/elastic_collisions.py new file mode 100644 index 0000000..e13140d --- /dev/null +++ b/build/lib/pycontest/elastic_collisions.py @@ -0,0 +1,48 @@ +import numpy as np + +def collision_1d(v1_i, v2_i, m1=1, m2=1): + """ 1d elastic collisions + + Arguments: + v1_1, v2_i: initial values of particles velocities + m1, m2: masses + + Returns: + final velocities after collisions + """ + + v1_f = (m1 - m2) / (m1 + m2) * v1_i + 2 * m2 / (m1 + m2) * v2_i + v2_f = (m2 - m1) / (m1 + m2) * v2_i + 2 * m1 / (m1 + m2) * v1_i + return v1_f, v2_f + + +def collision_2d(v1, v2, r12, m1, m2): + """ 2d elastic collisions + + Arguments: + v1_1, v2_i: initial values of particles velocities + m1, m2: masses + r12: vector between particles: r2-r1 + + Returns: + final velocities after collisions + """ + + if np.linalg.norm(r12) == 0: + raise Exception("two particles are exactly in the same place, " + "consider using smaller time step") + n_vect = r12 / np.linalg.norm(r12) + + # scalars + v1n = np.dot(v1, n_vect) + v2n = np.dot(v2, n_vect) + v1s = v1 - v1n * n_vect + v2s = v2 - v2n * n_vect + + v1n_f, v2n_f = collision_1d(v1n, v2n, m1, m2) + v1n_f *= n_vect + v2n_f *= n_vect + + v1_f = v1s + v1n_f + v2_f = v2s + v2n_f + return v1_f, v2_f \ No newline at end of file diff --git a/build/lib/pycontest/movies.py b/build/lib/pycontest/movies.py new file mode 100644 index 0000000..a3da8a5 --- /dev/null +++ b/build/lib/pycontest/movies.py @@ -0,0 +1,126 @@ +import matplotlib +matplotlib.use("TkAgg") +import matplotlib.pyplot as plt +import matplotlib.animation as anm +import numpy as np +import math as mt +import pdb + +""" +Matplotlib Animation Example + +author: Jake Vanderplas +email: vanderplas@astro.washington.edu +website: http://jakevdp.github.com +license: BSD +Please feel free to use and modify this, but keep the above information. Thanks! +""" +class Movie_2d: + + def __init__(self, fun, dt, t_max, loc, vel, domain, mass, radius=10): + self.dt = dt + self.fps = int(mt.floor(1 / dt)) + self.n_max = int(mt.floor(t_max / dt)) + self.loc = np.copy(loc) + self.vel = np.copy(vel) + self.domain = domain + self.mass = mass + self.radius = radius + self.fun = fun + + self.fig = plt.figure() + self.fig.subplots_adjust(left=0.1, bottom=0.1, right=.9, top=.9) + + self.ax = self.fig.add_subplot(111, aspect='equal', autoscale_on=False, xlim=(domain[0]), ylim=domain[1]) + for el in range(loc.shape[0]): + self.circ = plt.Circle(self.loc[el, :], self.radius, color='r', fill=False, linewidth=2) + self.ax.add_patch(self.circ) + + self.frame = plt.Rectangle([0, 0], domain[0][1] - domain[0][0], domain[1][1] - domain[1][0], fc='none') + self.ax.add_patch(self.frame) + + def generate_frame(self, i): + # pdb.set_trace() + loc, vel = self.fun(self.dt, self.mass, self.radius, self.loc, self.vel, self.domain) + #pdb.set_trace() + for el in range(loc.shape[0]): + self.circ.set_center(loc[el, :]) + return self.circ, self.frame + + def animate(self, name): + animation = anm.FuncAnimation(self.fig, self.generate_frame, frames=self.n_max, interval=10, blit=True) + + fps = self.fps + try: + animation.save(name + '.mp4', fps=fps, extra_args=['-vcodec', 'libx264']) + except ValueError: + print("\n mp4 file can not be created, will try to create html") + try: + animation.save(name + '.html', fps=fps, extra_args=['-vcodec', 'libx264']) + except: + print("\n html also can not be created") + +class Plot_2d: + + def __init__(self, dt, t_max, loc1_0, loc2_0, vel1_0, vel2_0, + domain, radius, m1=1, m2=1, nr=4, nc=5): + self.fig = plt.figure() + self.nr = nr + self.nc = nc + self.dt = dt + self.t_max = t_max + self.loc1_0 = loc1_0 + self.loc2_0 = loc2_0 + self.vel1_0 = vel1_0 + self.vel2_0 = vel2_0 + self.domain = domain + self.radius = radius + self.m1 = m1 + self.m2 = m2 + + + def simulation(self): + tt = 0 + ii = 0 + loc1, loc2 = self.loc1_0, self.loc2_0 + vel1, vel2 = self.vel1_0, self.vel2_0 + while tt <= t_max: + [loc1, loc2], [vel1, vel2] = \ + simulation_step(dt=self.dt, m1=self.m1, m2=self.m2, radius=self.radius, + loc1_0=loc1, loc2_0=loc2, v1_0=vel1, v2_0=vel2, + domain=self.domain) + tt += dt + if ii % 5 == 0: + #pdb.set_trace() + self.plotting(ii/5+1, loc1.copy(), loc2.copy(), vel1.copy(), vel2.copy()) + ii += 1 + plt.show() + + def plotting(self, ii, loc1, loc2, vel1, vel2): + ax = self.fig.add_subplot(self.nr, self.nc, ii, aspect='equal') + # ax = plt.gca() + # ax.quiver(r1[0], r1[1], r12[0], r12[1], angles='xy', scale_units='xy', scale=1, color="k", width=1.e-3) + # ax.quiver(0, 0, r1[0], r1[1], angles='xy', scale_units='xy', scale=1, color="k", width=1.e-3) + # ax.quiver(0, 0, r2[0], r2[1], angles='xy', scale_units='xy', scale=1, color="k", width=1.e-3) + # ax.quiver(r1[0], r1[1], v1[0], v1[1], angles='xy', scale_units='xy', scale=1, color="g") + # ax.quiver(r2[0], r2[1], v2[0], v2[1], angles='xy', scale_units='xy', scale=1, color="g") + # ax.quiver(r1[0], r1[1], v1n[0], v1n[1], angles='xy', scale_units='xy', scale=1, color="y") + # ax.quiver(r2[0], r2[1], v2n[0], v2n[1], angles='xy', scale_units='xy', scale=1, color="y") + # ax.quiver(r1[0], r1[1], v1s[0], v1s[1], angles='xy', scale_units='xy', scale=1, color="y") + # ax.quiver(r2[0], r2[1], v2s[0], v2s[1], angles='xy', scale_units='xy', scale=1, color="y") + # + # ax.quiver(r1[0], r1[1], u1[0], u1[1], angles='xy', scale_units='xy', scale=1, color="b") + # ax.quiver(r2[0], r2[1], u2[0], u2[1], angles='xy', scale_units='xy', scale=1, color="b") + # ax.quiver(r1[0], r1[1], u1n[0], u1n[1], angles='xy', scale_units='xy', scale=1, color="c") + # ax.quiver(r2[0], r2[1], u2n[0], u2n[1], angles='xy', scale_units='xy', scale=1, color="c") + # ax.quiver(r1[0], r1[1], u1s[0], u1s[1], angles='xy', scale_units='xy', scale=1, color="c") + # ax.quiver(r2[0], r2[1], u2s[0], u2s[1], angles='xy', scale_units='xy', scale=1, color="c") + + # ax.quiver(X, Y, V, U, angles='xy', scale_units='xy', scale=1, color="b") + c1 = plt.Circle(loc1, self.radius, color='r', fill=False) + c2 = plt.Circle(loc2, self.radius, color='r', fill=False) + ax.set_xlim([0, 100]) + ax.set_ylim([0, 100]) + # ax.set_aspec("equal") + ax.add_patch(c1) + ax.add_patch(c2) diff --git a/build/lib/pycontest/simulation.py b/build/lib/pycontest/simulation.py new file mode 100644 index 0000000..9ccd19a --- /dev/null +++ b/build/lib/pycontest/simulation.py @@ -0,0 +1,140 @@ +import numpy as np +import math as mt + +from scipy.spatial.distance import pdist, squareform + +from pycontest import elastic_collisions as ec +import pycontest.movies as mv +from pycontest.transport import transport + + +def simulation_step(dt, mass, radius, loc, vel, domain): + """ calculating positions and velocities of particles after one step + + Arguments: + dt: time steps + mass: masses + radius: radius for all balls + loc: locations of balls + vel: 2d velocities of balls + domain: domain for the box + Returns: + positions and velocities of particles after one step + """ + + + # transport + loc = transport(loc, vel, dt) + + r12_x = squareform(pdist(loc, lambda r1, r2: r1[0] - r2[0])) + r12_y = squareform(pdist(loc, lambda r1, r2: r1[1] - r2[1])) + r12 = np.stack((r12_x, r12_y), axis=-1) + dist = np.apply_along_axis(np.linalg.norm, 2, r12) + + # finding particles that are colliding + ind1, ind2 = np.where(dist < 2 * radius) + unique = (ind1 < ind2) + + ind1 = ind1[unique] + ind2 = ind2[unique] + + # collisions + for id1, id2 in zip(ind1, ind2): + vel[id1], vel[id2] = ec.collision_2d(vel[id1], vel[id2], r12[id1, id2], mass[id1], mass[id2]) + + # find outside domain points ... + out_left = (loc[:, 0] < domain[0][0] + radius) + out_right = (loc[:, 0] > domain[0][1] - radius) + out_top = (loc[:, 1] < domain[1][0] + radius) + out_bottom = (loc[:, 1] > domain[1][1] - radius) + + # ... change their velocities + vel[out_left | out_right, 0] *= -1 + vel[out_top | out_bottom, 1] *= -1 + + return loc, vel + + +def simulation(t_max, dt, mass, radius, loc_0, vel_0, domain, t0=0): + """ calculating positions and velocities of particles after t_max + + Arguments: + t_max: length of the simulation + dt: time steps + mass: masses + radius: radius for all balls + loc_0: initial locations of balls + vel_0: 2d initial velocities of balls + t0: initial time (0 as default) + domain: domain for the box + Returns: + positions and velocities of particles at t_max + """ + t = t0 + loc = np.copy(loc_0) + vel = np.copy(vel_0) + while(t