#!/usr/bin/env python
# coding: utf-8

# In[1]:


get_ipython().run_line_magic('pylab', 'inline')


# # Efusión

# In[2]:


#Load the simulation
from pythonmd.simcore import MDSimulationSession


# In[3]:


get_ipython().run_cell_magic('file', 'effusion', '[MD Simulation]\nnumber_of_molecules: 400\nbox_width: 600\nbox_height: 300\nwall: effusion\nwall_pos: 300\nwall_hole: 0\n')


# In[6]:


time = []
n_left = []
n_right = []
t_left = []
t_right = []

session = MDSimulationSession('effusion')
session.run_until_t(4)
session.run(100)
session.run_until_t(2)
session.run_until_equilibrium(delta_t=1e-05)
current_time, t_l, _, _, n_l = session.get_left_stats()
_, t_r, _, _, n_r = session.get_right_stats()

time.append(current_time)
n_left.append(n_l)
n_right.append(n_r)
t_left.append(t_l)
t_right.append(t_r)

session.change_hole(10)
session.reset_statistics()

for _ in range(60000):
    session.run(1)
    current_time, t_l, _, _, n_l = session.get_left_stats()
    _, t_r, _, _, n_r = session.get_right_stats()
    
    time.append(current_time)
    n_left.append(n_l)
    n_right.append(n_r)
    t_left.append(t_l)
    t_right.append(t_r)


plt.plot(time, n_left, label='left')
plt.plot(time, n_right, color='red', label='right')
plt.xlabel("Time ($\sigma\sqrt{m/\epsilon}$)")
plt.ylabel("Number of Molecules")
plt.legend(loc='best')


# In[7]:


plt.plot(time, t_left, label='left')
plt.plot(time, t_right, color='red', label='right')
plt.xlabel("Time ($\sigma\sqrt{m/\epsilon}$)")
plt.ylabel("Temperature ($\epsilon/k_b$)")
plt.legend(loc='best')


# In[27]:


average_n_left = (sum(n_left[32000:])/len(n_left[32000:]))
average_t_left = (sum(t_left[32000:])/len(t_left[32000:]))
average_n_right = (sum(n_right[32000:])/len(n_right[32000:]))
average_t_right = (sum(t_right[32000:])/len(t_right[32000:]))

print("Initial time of equilibrium: %i" % time[32000])
print("Average n left: %.4f" % average_n_left)
print("Average T left: %.4f" % average_t_left)
print("Average n right: %.4f" % average_n_right)
print("Average T right: %.4f" % average_t_right)
print(average_n_left*math.sqrt(average_t_left)/(average_n_right*math.sqrt(average_t_right)))