<reponame>alexeltsar/Gas-Simulation<filename>Gas_classes.py
import numpy as np
import math, time
import random
import matplotlib.pyplot as plt
from numpy import linalg as LA
import mpl_toolkits.mplot3d.axes3d as p3
from matplotlib.animation import FuncAnimation
import matplotlib.patches as mpatches
from scipy import stats

import os
import sqlite3

k_B = 1.380648e-23 #Boltzmann constant (J/K)

current_milli_time = lambda: int(round(time.time() * 1000))

dir_principal = os.getcwd()
carpeta_data = dir_principal + '\Data'

###########################################################
#														  #
#						FUNCTIONS						  #
#													  	  #
###########################################################

def mod(v):
    """
        Computes the squared sum over the last axis of the numpy.ndarray v
    """
    return np.sum(v * v, axis=-1)

def calculate_T(v, m, N):
	suma = np.sum(v * v)

	temperature = suma * m / (3*N * k_B)

	return temperature

def calculate_P(N, x, T, dt, p):
	average_p = np.sum(p) / len(p)
	Pressure = average_p / (dt * 6 * x**2)

	theory = N * k_B * T / (x**3)

	return (Pressure, theory)

def vel_hist(x, v, v_max, dv, m, T):

	def fmb(m, T, v):
		return (m/(2.*np.pi*k_B*T))**1.5 * 4.*np.pi*v*v * np.exp((-m*v*v)/(2.*k_B*T))  #Maxwell-Boltzman distribution
 
	v = np.arange(0.0, v_max, dv)

	plt.plot(v, fmb(m, T, v), 'r-')
	n, bins, patches = plt.hist(x, 40, normed=1, facecolor='g', alpha=0.5)
	plt.xlabel('Speed -- v (m/s)')
	plt.ylabel('Probability -- f(v)')
	plt.title('Histogram of Maxwell-Boltzmann speed distribution')
	plt.grid(True)
	
	os.chdir(carpeta_data)
	plt.savefig('Velocity_hist.png')
	plt.gcf().clear()

def create_database():
	database = sqlite3.connect("gas_database.db")
	cursor = database.cursor()

	tablas = [
			"""
				CREATE TABLE IF NOT EXISTS ideal_gas(
					temperature REAL NOT NULL,
					theo_pressure REAL NOT NULL,
					exp_pressure REAL NOT NULL,
					rel_error REAL NOT NULL
				);
			"""
		]
	for tabla in tablas:
		cursor.execute(tabla);

	tablas = [
			"""
				CREATE TABLE IF NOT EXISTS real_gas(
					temperature REAL NOT NULL,
					theo_pressure REAL NOT NULL,
					exp_pressure REAL NOT NULL,
					rel_error_pressure REAL NOT NULL,
					theo_avg_vel REAL NOT NULL,
					exp_avg_vel REAL NOT NULL,
					rel_error_vel REAL NOT NULL
				);
			"""
		]
	for tabla in tablas:
		cursor.execute(tabla);

	database.commit()

def fill_database(name, table, interrogants, values):
	database = sqlite3.connect('%s.db' % name)
	cursor = database.cursor()

	sentencia = "INSERT INTO %s VALUES %s" % (table, interrogants) 
	cursor.execute(sentencia, values)
	database.commit()
	
def read_database(name, table, name_2, order):
	database = sqlite3.connect('%s.db' % name)
	cursor = database.cursor()
	sentencia = "SELECT * FROM '%s' ORDER BY %s %s;" % (table, name_2, order)
 
	cursor.execute(sentencia)
	
	lines = cursor.fetchall()

	cursor.close()
	
	return lines

###########################################################
#														  #
#						CLASSES						  	  #
#													  	  #
###########################################################

class IdealGas(object):
	def __init__(self, N, L, T, m, R, dt, T_MAX):
		self.N = N #Number of particles
		self.L = L #Length of the box
		self.T = T #Temperature of the system
		self.m = m #Particle mass (Oxygen  2.66e-23) (He 6.64e-24)
		self.R = R #Particle radius
		self.dt = dt #Step size
		self.T_MAX = T_MAX #Maximum simulation time

		self.Nt = int(self.T_MAX / self.dt) #Number of simulation steps/iterations

		self.r = np.zeros((self.N, 3)) #Array of vector positions
		self.v = np.zeros((self.N, 3)) #Array of vector velocities
		self.p = np.zeros(int(self.Nt)) # Array of momentum exchanged with the wall
		self.CM = (0,0,0) #Center of masses

		self.CM_t = []
		
	def init_particles(self):

		self.r = np.random.rand(self.N, 3) * 2 * (self.L/2 - self.R) - (self.L/2 - self.R)

		#vmed = np.sqrt(8 * k_B * T / (np.pi * m))
		vrms = np.sqrt(3*k_B * self.T / self.m)

		v_polar = np.random.random((self.N, 2))

		self.v[:,0] = (np.sin(v_polar[:,0] * np.pi) * np.cos(v_polar[:,1] * 2 * np.pi)) 
		self.v[:,1] = (np.sin(v_polar[:,0] * np.pi) * np.sin(v_polar[:,1] * 2 * np.pi)) 
		self.v[:,2] = np.cos(v_polar[:,0] * np.pi)

		self.v *= vrms

	def check_walls(self, i, x): #Checks for particle-wall collisions and return the exchanged momentum
		pex = 0
		for k in range(3):
			if abs(self.r[i][k]) + self.R > x/2:
				self.v[i][k] *= -1
				if self.r[i][k] < 0:
					self.r[i][k] = -x/2 + self.R
				else :
					self.r[i][k] = x/2 - self.R

				pex += 2 * self.m * abs(self.v[i][k])
		return pex

	def simulation(self):
		self.init_particles()

		for k in range(self.Nt):
			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]

			self.CM_t.append(np.sum(self.r, axis=0) / self.N)

		return self.p

	def isothermal_expansion(self, pull_period, L0, dL):
		self.init_particles()

		Vs = []
		Vs_discrete = []
		theo_pr_discrete = []
		time = []
		Exp_pr = []
		Theo_pr = []


		Vs_discrete.append(L0**3)	
		theo_pr_discrete.append(self.N * k_B * self.T / (L0**3))

		self.L = L0

		for k in range(self.Nt):
			if k % (pull_period/self.dt) == 0 and k > 0:
				self.L = self.L + dL
				Vs_discrete.append(self.L**3)
				theo_pr_discrete.append(self.N * k_B* self.T/(self.L**3))

			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]

			if k > 1:

				Exp_P = self.p[k] / (self.dt * 6 * self.L**2)
				Theo_P = self.N * k_B * self.T / (self.L**3)

				Exp_pr.append(Exp_P)
				Theo_pr.append(Theo_P)
				Vs.append(self.L**3)
				time.append(k*self.dt)

		if pull_period > self.dt :		
			a = int(self.T_MAX/pull_period) #Times the pressure will remain constant (as V is doesn't change)

			average_t = [] #average of the pressure in this times

			for h in range(a):
				avrg = 0
				realtime = int(pull_period/self.dt)
				for s in range(h*realtime, (h+1)*realtime-2):
					avrg = avrg + Exp_pr[s]
				average_t.append(avrg/realtime)

		plt.plot(Vs, Theo_pr, ':o', label = 'Theoretical pressure')
		plt.plot(Vs, Exp_pr, 'x', label = 'Experimental pressure')

		if pull_period > self.dt :
			plt.plot(Vs_discrete, average_t, ls = 'none', marker = 'o', label = 'Average experimental pressure')

		plt.title('Isothermal compression (T=%.2f K)' % self.T)
		plt.xlabel('Volume ($m^3$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()
		
		os.chdir(carpeta_data)
		plt.savefig('Expansion_volume.png')
		plt.gcf().clear()

		plt.plot(time, Theo_pr, '-', label = 'Theoretical pressure')
		plt.plot(time, Exp_pr, 'x', label = 'Experimental pressure')

		if pull_period > self.dt :
			plt.plot(np.linspace(0+ pull_period/2, self.T_MAX - pull_period/2, (a)), average_t, ls = 'none', marker='o', label = 'Average experimental pressure')

		plt.title('Isothermal compression (T=%.2f K)' % self.T)
		plt.xlabel('time ($s$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()
		
		os.chdir(carpeta_data)
		plt.savefig('Expansion_time.png')
		plt.gcf().clear()			

	def isochoric_heating(self, heat_period, T0, dT):	
		T_discrete = []
		T_discrete.append(T0)
		counter =0

		Exp_pr = []
		Theo_pr = []
		temperatures = []
		time = []

		self.init_particles()

		for k in range(self.Nt):

			if k % (heat_period / self.dt) == 0 and k > 0:
				self.T = self.T + dT
				T_discrete.append(self.T)

				vmed = np.sqrt(8 * k_B * self.T / (np.pi * self.m)) #Re-init particles velocities for each T
				v_polar = np.random.random((self.N, 2))

				self.v[:,0] = (np.sin(v_polar[:,0] * np.pi) * np.cos(v_polar[:,1] * 2 * np.pi)) 
				self.v[:,1] = (np.sin(v_polar[:,0] * np.pi) * np.sin(v_polar[:,1] * 2 * np.pi)) 
				self.v[:,2] = np.cos(v_polar[:,0] * np.pi)

				self.v *= vmed

			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]

			Exp_P = self.p[k]/(self.dt * 6 * self.L**2)
			Theo_P = self.N * k_B * self.T / (self.L**3)

			Exp_pr.append(Exp_P)
			Theo_pr.append(Theo_P)
			temperatures.append(self.T)
			time.append(k * self.dt)

		if heat_period > self.dt:
			a = int(self.T_MAX / heat_period) #Times the pressure will remain constant (as T doesn't change)

			average_t = [] #average of the pressure in this times

			for h in range(a):
				avrg = 0
				realtime = int(heat_period / self.dt)
				for s in range(h*realtime, (h+1)*realtime-2):
					avrg = avrg + Exp_pr[s]
				average_t.append(avrg/realtime)

		plt.plot(temperatures, Theo_pr, ':.', label = 'Theoretical pressure')
		plt.plot(temperatures, Exp_pr, 'x', label = 'Experimental pressure')

		if heat_period > self.dt:
			plt.plot(T_discrete, average_t, '-o', label = 'Average experimental pressure')

		plt.title('Isochoric process (V=%.2f $m^3$)'%(self.L**3))
		plt.xlabel('Temperature ($K$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()

		os.chdir(carpeta_data)
		plt.savefig('Heating_temperature.png')
		plt.gcf().clear()

		plt.plot(time, Theo_pr, '-', label = 'Theoretical pressure')
		plt.plot(time, Exp_pr, 'x', label = 'Experimental pressure')

		if heat_period > self.dt:
			plt.plot(np.linspace(0+ heat_period/2, self.T_MAX - heat_period/2, (a)), average_t, ls = 'none', marker='o', label = 'Average experimental pressure')

		plt.title('Isochoric process (V=%.2f $m^3$)'%(self.L**3))
		plt.xlabel('time ($s$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()

		os.chdir(carpeta_data)
		plt.savefig('Heating_time.png')
		plt.gcf().clear()

	def animation(self):

		self.init_particles()

		def update(t, lines):
			k = int(t / self.dt)

			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]

			self.CM = np.sum(self.r, axis=0) / self.N
			self.CM_t.append(self.CM)

			lines[0].set_data(self.r[:,0], self.r[:,1])
			lines[0].set_3d_properties(self.r[:,2])

			lines[1].set_data([self.CM[0]], [self.CM[1]])
			lines[1].set_3d_properties([self.CM[2]])

			return lines

		# Attaching 3D axis to the figure
		fig = plt.figure()
		ax = p3.Axes3D(fig)


		# Setting the axes properties
		ax.set_xlim3d([-self.L/2, self.L/2])
		ax.set_xlabel('X')

		ax.set_ylim3d([-self.L/2, self.L/2])
		ax.set_ylabel('Y')

		ax.set_zlim3d([-self.L/2, self.L/2])
		ax.set_zlabel('Z')

		T1 = calculate_T(self.v, self.m, self.N)

		lines = []
		lines.append(ax.plot(self.r[:,0], self.r[:,1], self.r[:,2], ls='None', marker='.', label = 'Particles at $T_0=%.2f\,K$'%T1)[0])
		lines.append(ax.plot([self.CM[0]], [self.CM[1]], [self.CM[2]], marker='o', color='r', label='Center of masses')[0])

		ani = FuncAnimation(fig, update, fargs=(lines,), frames=np.linspace(0, self.T_MAX-self.dt, int(self.T_MAX / self.dt)),
		                    blit=True, interval=1, repeat=False)
		#ani.save('animation.mp4', fps=20, writer="ffmpeg", codec="libx264")
		plt.legend(loc = 'upper left')
		plt.show()

	def experiments(self):

		Exp_P, Theo_P = calculate_P(self.N, self.L, self.T, self.dt, self.p)
		T1 = calculate_T(self.v, self.m, self.N)

		A = 6*self.L**2

		theo_P_2 = self.N * k_B * T1 / (self.L**3)

		rel_error = abs(Theo_P - Exp_P) / Theo_P * 100
		rel_error2 = abs(theo_P_2 - Exp_P) / theo_P_2 * 100

		os.chdir(carpeta_data)
		fill_database('gas_database', 'ideal_gas', '(?,?,?,?)', [T1, Theo_P, Exp_P, rel_error])

		#Pressure plot
		theo_array = []
		exp_array = []

		for i in range(self.Nt):
			theo_array.append(Theo_P)
			exp_array.append(Exp_P)

		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), self.p / (self.dt * A), ':', label = 'Experimental pressure')
		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), theo_array, label = 'Theoretical pressure')
		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), exp_array, label = 'Average experimental pressure')
		plt.title('Pressure of the system')
		plt.xlabel('t(s)')
		plt.ylabel('P(Pa)')
		plt.legend()

		os.chdir(carpeta_data)
		plt.savefig('Pressure.png')
		plt.gcf().clear()

		#CM PLOTS
		transposed = np.transpose(self.CM_t)
		t = np.linspace(0, self.T_MAX - self.dt, len(self.CM_t))

		std_x = np.std(transposed[:][0])
		std_y = np.std(transposed[:][1])
		std_z = np.std(transposed[:][2])

		mean_x = np.mean(transposed[:][0])
		mean_y = np.mean(transposed[:][1])
		mean_z = np.mean(transposed[:][2])

		plt.subplot(2, 2, 1)
		plt.plot(t, transposed[:][0]) #CM_x
		meanx = plt.plot(t, np.linspace(mean_x, mean_x, len(self.CM_t)), color='red', label='Mean') #mean
		sdx = plt.plot(t, np.linspace(std_x, std_x, len(self.CM_t)), ls = ':', color='green', label=r'$\pm$ SD') #+ standard deviation
		plt.plot(t, np.linspace(-std_x, -std_x, len(self.CM_t)), ls = ':', color='green') #- standard deviation
		plt.xlabel('t')
		plt.ylabel('x')

		plt.subplot(2, 2, 2)
		plt.plot(t, transposed[:][1]) #CM_y
		plt.plot(t, np.linspace(mean_y, mean_y, len(self.CM_t)), color='red', label='Mean') #mean
		plt.plot(t, np.linspace(std_y, std_y, len(self.CM_t)), ls = ':', color='green', label=r'$\pm$ SD') #+ standard deviation
		plt.plot(t, np.linspace(-std_y, -std_y, len(self.CM_t)), ls = ':', color='green') #- standard deviation
		plt.xlabel('t')
		plt.ylabel('y')

		plt.subplot(2, 2, 3)
		plt.plot(t, transposed[:][2]) #CM_z
		plt.plot(t, np.linspace(mean_z, mean_z, len(self.CM_t)), color='red', label='Mean') #mean
		plt.plot(t, np.linspace(std_z, std_z, len(self.CM_t)), ls = ':', color='green', label=r'$\pm$ SD') #+ standard deviation
		plt.plot(t, np.linspace(-std_z, -std_z, len(self.CM_t)), ls = ':', color='green') #- standard deviation
		plt.xlabel('t')
		plt.ylabel('z')

		plt.suptitle(r'$\vec{r}_{CM}$ vs t')
		plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
		
		os.chdir(carpeta_data)
		plt.savefig('CM_plots.png')
		plt.gcf().clear()

class RealGas(object):
	def __init__(self, N, L, T, m, R, dt, T_MAX):
		self.N = N #Number of particles
		self.L = L #Length of the box
		self.T = T #Temperature of the system
		self.m = m #Particle mass (Oxygen  2.66e-23) (He 6.64e-24)
		self.R = R #Particle radius
		self.dt = dt #Step size
		self.T_MAX = T_MAX #Maximum simulation time

		self.Nt = int(self.T_MAX / self.dt) #Number of simulation steps/iterations

		self.r = np.zeros((self.N, 3)) #Array of vector positions
		self.v = np.zeros((self.N, 3)) #Array of vector velocities
		self.p = np.zeros(int(self.Nt)) # Array of momentum exchanged with the wall
		self.CM = (0,0,0) #Center of masses

		self.hist = []
		self.CM_t = []

	def init_particles(self):

		self.r = np.random.rand(self.N, 3) * 2 * (self.L/2 - self.R) - (self.L/2 - self.R)

		#vmed = np.sqrt(8 * k_B * T / (np.pi * m))
		vrms = np.sqrt(3*k_B * self.T / self.m)

		v_polar = np.random.random((self.N, 2))

		self.v[:,0] = (np.sin(v_polar[:,0] * np.pi) * np.cos(v_polar[:,1] * 2 * np.pi)) 
		self.v[:,1] = (np.sin(v_polar[:,0] * np.pi) * np.sin(v_polar[:,1] * 2 * np.pi)) 
		self.v[:,2] = np.cos(v_polar[:,0] * np.pi)

		self.v *= vrms

	def check_walls(self, i, x): #Checks for particle-wall collisions and return the exchanged momentum
		pex = 0
		for k in range(3):
			if abs(self.r[i][k]) + self.R > x/2:
				self.v[i][k] *= -1
				if self.r[i][k] < 0:
					self.r[i][k] = -x/2 + self.R
				else :
					self.r[i][k] = x/2 - self.R

				pex += 2 * self.m * abs(self.v[i][k])
		return pex

	def check_collision(self):  #Check for collisions with other particles of te same gas

		dists = np.sqrt(mod(self.r - self.r[:,np.newaxis])) 
		cols2 = (0 < dists) & (dists < 2*self.R)
		idx_i, idx_j = np.nonzero(cols2)
		# ***possibility to simplify this *** #
		for i, j in zip(idx_i, idx_j):
		    if j < i:
		        # skip duplications and same particle
		        continue 

		    rij = self.r[i] - self.r[j]
		    d = mod(rij)
		    vij = self.v[i] - self.v[j]
		    dv = np.dot(vij, rij) * rij / d
		    self.v[i] -= dv
		    self.v[j] += dv
			
		    # update the positions so they are no longer in contact
		    self.r[i] += self.dt * self.v[i]
		    self.r[j] += self.dt * self.v[j]

	def simulation(self):
		self.init_particles()

		for k in range(self.Nt):
			self.check_collision()

			modv = 0
			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]
				modv = modv + LA.norm(self.v[i])

			self.CM_t.append(np.sum(self.r, axis=0) / self.N)

			average_v = modv / self.N
			self.hist.append(average_v)

		return self.p

	def isothermal_expansion(self, pull_period, L0, dL):
		self.init_particles()

		Vs = []
		Vs_discrete = []
		theo_pr_discrete = []
		time = []
		Exp_pr = []
		Theo_pr = []


		Vs_discrete.append(L0**3)	
		theo_pr_discrete.append(self.N * k_B * self.T / (L0**3))

		self.L = L0

		for k in range(self.Nt):
			self.check_collision()

			if k % (pull_period/self.dt) == 0 and k > 0:
				self.L = self.L + dL
				Vs_discrete.append(self.L**3)
				theo_pr_discrete.append(self.N * k_B* self.T/(self.L**3))

			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]

			if k > 1:

				Exp_P = self.p[k] / (self.dt * 6 * self.L**2)
				Theo_P = self.N * k_B * self.T / (self.L**3)

				Exp_pr.append(Exp_P)
				Theo_pr.append(Theo_P)
				Vs.append(self.L**3)
				time.append(k*self.dt)

		if pull_period > self.dt :		
			a = int(self.T_MAX/pull_period) #Times the pressure will remain constant (as V is doesn't change)

			average_t = [] #average of the pressure in this times

			for h in range(a):
				avrg = 0
				realtime = int(pull_period/self.dt)
				for s in range(h*realtime, (h+1)*realtime-2):
					avrg = avrg + Exp_pr[s]
				average_t.append(avrg/realtime)

		plt.plot(Vs, Theo_pr, ':o', label = 'Theoretical pressure')
		plt.plot(Vs, Exp_pr, 'x', label = 'Experimental pressure')

		if pull_period > self.dt :
			plt.plot(Vs_discrete, average_t, ls = 'none', marker = 'o', label = 'Average experimental pressure')

		plt.title('Isothermal compression (T=%.2f K)' % self.T)
		plt.xlabel('Volume ($m^3$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()
		
		os.chdir(carpeta_data)
		plt.savefig('Expansion_volume.png')
		plt.gcf().clear()

		plt.plot(time, Theo_pr, '-', label = 'Theoretical pressure')
		plt.plot(time, Exp_pr, 'x', label = 'Experimental pressure')

		if pull_period > self.dt :
			plt.plot(np.linspace(0+ pull_period/2, self.T_MAX - pull_period/2, (a)), average_t, ls = 'none', marker='o', label = 'Average experimental pressure')

		plt.title('Isothermal compression (T=%.2f K)' % self.T)
		plt.xlabel('time ($s$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()
		
		os.chdir(carpeta_data)
		plt.savefig('Expansion_time.png')
		plt.gcf().clear()			

	def isochoric_heating(self, heat_period, T0, dT):	
		T_discrete = []
		T_discrete.append(T0)
		counter =0

		Exp_pr = []
		Theo_pr = []
		temperatures = []
		time = []

		self.init_particles()

		for k in range(self.Nt):
			self.check_collision()

			if k % (heat_period / self.dt) == 0 and k > 0:
				self.T = self.T + dT
				T_discrete.append(self.T)

				vmed = np.sqrt(8 * k_B * self.T / (np.pi * self.m)) #Re-init particles velocities for each T
				v_polar = np.random.random((self.N, 2))

				self.v[:,0] = (np.sin(v_polar[:,0] * np.pi) * np.cos(v_polar[:,1] * 2 * np.pi)) 
				self.v[:,1] = (np.sin(v_polar[:,0] * np.pi) * np.sin(v_polar[:,1] * 2 * np.pi)) 
				self.v[:,2] = np.cos(v_polar[:,0] * np.pi)

				self.v *= vmed

			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]

			Exp_P = self.p[k]/(self.dt * 6 * self.L**2)
			Theo_P = self.N * k_B * self.T / (self.L**3)

			Exp_pr.append(Exp_P)
			Theo_pr.append(Theo_P)
			temperatures.append(self.T)
			time.append(k * self.dt)

		if heat_period > self.dt:
			a = int(self.T_MAX / heat_period) #Times the pressure will remain constant (as T doesn't change)

			average_t = [] #average of the pressure in this times

			for h in range(a):
				avrg = 0
				realtime = int(heat_period / self.dt)
				for s in range(h*realtime, (h+1)*realtime-2):
					avrg = avrg + Exp_pr[s]
				average_t.append(avrg/realtime)

		plt.plot(temperatures, Theo_pr, ':.', label = 'Theoretical pressure')
		plt.plot(temperatures, Exp_pr, 'x', label = 'Experimental pressure')

		if heat_period > self.dt:
			plt.plot(T_discrete, average_t, '-o', label = 'Average experimental pressure')

		plt.title('Isochoric process (V=%.2f $m^3$)'%(self.L**3))
		plt.xlabel('Temperature ($K$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()

		os.chdir(carpeta_data)
		plt.savefig('Heating_temperature.png')
		plt.gcf().clear()

		plt.plot(time, Theo_pr, '-', label = 'Theoretical pressure')
		plt.plot(time, Exp_pr, 'x', label = 'Experimental pressure')

		if heat_period > self.dt:
			plt.plot(np.linspace(0+ heat_period/2, self.T_MAX - heat_period/2, (a)), average_t, ls = 'none', marker='o', label = 'Average experimental pressure')

		plt.title('Isochoric process (V=%.2f $m^3$)'%(self.L**3))
		plt.xlabel('time ($s$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()

		os.chdir(carpeta_data)
		plt.savefig('Heating_time.png')
		plt.gcf().clear()


	def animation(self):

		self.init_particles()

		def update(t, lines):
			k = int(t / self.dt)
			self.check_collision(self.r, self.v)

			modv = 0
			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]
				modv = modv + LA.norm(self.v[i])

			self.CM = np.sum(self.r, axis=0) / self.N
			self.CM_t.append(self.CM)

			average_v = modv / self.N
			self.hist.append(average_v)

			lines[0].set_data(self.r[:,0], self.r[:,1])
			lines[0].set_3d_properties(self.r[:,2])

			lines[1].set_data([self.CM[0]], [self.CM[1]])
			lines[1].set_3d_properties([self.CM[2]])

			return lines

		# Attaching 3D axis to the figure
		fig = plt.figure()
		ax = p3.Axes3D(fig)


		# Setting the axes properties
		ax.set_xlim3d([-self.L/2, self.L/2])
		ax.set_xlabel('X')

		ax.set_ylim3d([-self.L/2, self.L/2])
		ax.set_ylabel('Y')

		ax.set_zlim3d([-self.L/2, self.L/2])
		ax.set_zlabel('Z')

		T1 = calculate_T(self.v, self.m, self.N)

		lines = []
		lines.append(ax.plot(self.r[:,0], self.r[:,1], self.r[:,2], ls='None', marker='.', label = 'Particles at $T_0=%.2f\,K$'%T1)[0])
		lines.append(ax.plot([self.CM[0]], [self.CM[1]], [self.CM[2]], marker='o', color='r', label='Center of masses')[0])

		ani = FuncAnimation(fig, update, fargs=(lines,), frames=np.linspace(0, self.T_MAX-self.dt, self.Nt),
		                    blit=True, interval=1, repeat=False)
		#ani.save('animation.mp4', fps=20, writer="ffmpeg", codec="libx264")
		plt.legend(loc = 'upper left')
		plt.show()

	def experiments(self):

		Exp_P, Theo_P = calculate_P(self.N, self.L, self.T, self.dt, self.p)
		T1 = calculate_T(self.v, self.m, self.N)

		A = 6*self.L**2

		theo_P_2 = self.N * k_B * T1 / (self.L**3)

		rel_error = abs(Theo_P - Exp_P) / Theo_P * 100
		rel_error2 = abs(theo_P_2 - Exp_P) / theo_P_2 * 100

		#PRESSURE PLOT
		theo_array = []
		exp_array = []

		for i in range(self.Nt):
			theo_array.append(Theo_P)
			exp_array.append(Exp_P)

		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), self.p / (self.dt * A), ':', label = 'Experimental pressure')
		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), theo_array, label = 'Theoretical pressure')
		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), exp_array, label = 'Average experimental pressure')
		plt.title('Pressure of the system')
		plt.xlabel('t(s)')
		plt.ylabel('P(Pa)')
		plt.legend()
		
		os.chdir(carpeta_data)
		plt.savefig('Pressure.png')
		plt.gcf().clear()

		#VELOCITY PLOT
		theo_array = []
		theo_array_rms = []

		vmed = np.sqrt(8 * k_B * self.T / (np.pi * self.m))
		vrms = np.sqrt(3 * k_B * self.T / self.m)

		exp_avg_vel = np.sum(self.hist) / len(self.hist)

		rel_error_vel = abs(vmed-exp_avg_vel)/vmed * 100

		fill_database('gas_database', 'real_gas', '(?,?,?,?,?,?,?)', [T1, Theo_P, Exp_P, rel_error, vmed, exp_avg_vel, rel_error_vel])

		plt.plot(np.linspace(0, self.T_MAX - self.dt, len(self.hist)), self.hist, '-', label = 'Experimental data')
		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), np.linspace(vmed, vmed, self.Nt), label = r'Maxwell-Boltzman $v_{med}$')
		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), np.linspace(vrms, vrms, self.Nt), label = r'Maxwell-Boltzman $v_{rms}$')
		plt.title('Average velocity ($T=%.2f\,K$)' % T1)
		plt.xlabel('t(s)')
		plt.ylabel('$v_{med}$(m/s)')
		plt.legend()
		
		os.chdir(carpeta_data)
		plt.savefig('Velocity.png')
		plt.gcf().clear()

		#VELOCITY HISTOGRAM

		modv = []
		for i in range(len(self.v)):
			mod = self.v[i][0]**2 + self.v[i][1]**2 + self.v[i][2]**2
			modv.append(np.sqrt(mod))

		v_max = 3*np.sqrt((3*k_B* self.T)/ self.m)
		dv = 0.01

		vel_hist(modv, self.v, v_max, dv, self.m, self.T)

		#CENTER OF MASSES MOVEMENT PLOTS

		transposed = np.transpose(self.CM_t)
		t = np.linspace(0, self.T_MAX - self.dt, len(self.CM_t))

		std_x = np.std(transposed[:][0])
		std_y = np.std(transposed[:][1])
		std_z = np.std(transposed[:][2])

		mean_x = np.mean(transposed[:][0])
		mean_y = np.mean(transposed[:][1])
		mean_z = np.mean(transposed[:][2])

		plt.subplot(2, 2, 1)
		plt.plot(t, transposed[:][0]) #CM_x
		meanx = plt.plot(t, np.linspace(mean_x, mean_x, len(self.CM_t)), color='red', label='Mean') #mean
		sdx = plt.plot(t, np.linspace(std_x, std_x, len(self.CM_t)), ls = ':', color='green', label=r'$\pm$ SD') #+ standard deviation
		plt.plot(t, np.linspace(-std_x, -std_x, len(self.CM_t)), ls = ':', color='green') #- standard deviation
		plt.xlabel('t')
		plt.ylabel('x')

		plt.subplot(2, 2, 2)
		plt.plot(t, transposed[:][1]) #CM_y
		plt.plot(t, np.linspace(mean_y, mean_y, len(self.CM_t)), color='red', label='Mean') #mean
		plt.plot(t, np.linspace(std_y, std_y, len(self.CM_t)), ls = ':', color='green', label=r'$\pm$ SD') #+ standard deviation
		plt.plot(t, np.linspace(-std_y, -std_y, len(self.CM_t)), ls = ':', color='green') #- standard deviation
		plt.xlabel('t')
		plt.ylabel('y')

		plt.subplot(2, 2, 3)
		plt.plot(t, transposed[:][2]) #CM_z
		plt.plot(t, np.linspace(mean_z, mean_z, len(self.CM_t)), color='red', label='Mean') #mean
		plt.plot(t, np.linspace(std_z, std_z, len(self.CM_t)), ls = ':', color='green', label=r'$\pm$ SD') #+ standard deviation
		plt.plot(t, np.linspace(-std_z, -std_z, len(self.CM_t)), ls = ':', color='green') #- standard deviation
		plt.xlabel('t')
		plt.ylabel('z')

		plt.suptitle(r'$\vec{r}_{CM}$ vs t')
		plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')

		os.chdir(carpeta_data)
		plt.savefig('CM_plots.png')
		plt.gcf().clear()


#Fer isothermal and isochoric processes a la classe de la interfase per utilitzar la barra de carrega

import numpy as np
import math, time
import random
import matplotlib.pyplot as plt
from numpy import linalg as LA
import mpl_toolkits.mplot3d.axes3d as p3
from matplotlib.animation import FuncAnimation
import matplotlib.patches as mpatches
from scipy import stats

import os
import sqlite3

k_B = 1.380648e-23 #Boltzmann constant (J/K)

current_milli_time = lambda: int(round(time.time() * 1000))

dir_principal = os.getcwd()
carpeta_data = dir_principal + '\Data'

###########################################################
#														  #
#						FUNCTIONS						  #
#													  	  #
###########################################################

def mod(v):
    """
        Computes the squared sum over the last axis of the numpy.ndarray v
    """
    return np.sum(v * v, axis=-1)

def calculate_T(v, m, N):
	suma = np.sum(v * v)

	temperature = suma * m / (3*N * k_B)

	return temperature

def calculate_P(N, x, T, dt, p):
	average_p = np.sum(p) / len(p)
	Pressure = average_p / (dt * 6 * x**2)

	theory = N * k_B * T / (x**3)

	return (Pressure, theory)

def vel_hist(x, v, v_max, dv, m, T):

	def fmb(m, T, v):
		return (m/(2.*np.pi*k_B*T))**1.5 * 4.*np.pi*v*v * np.exp((-m*v*v)/(2.*k_B*T))  #Maxwell-Boltzman distribution
 
	v = np.arange(0.0, v_max, dv)

	plt.plot(v, fmb(m, T, v), 'r-')
	n, bins, patches = plt.hist(x, 40, normed=1, facecolor='g', alpha=0.5)
	plt.xlabel('Speed -- v (m/s)')
	plt.ylabel('Probability -- f(v)')
	plt.title('Histogram of Maxwell-Boltzmann speed distribution')
	plt.grid(True)
	
	os.chdir(carpeta_data)
	plt.savefig('Velocity_hist.png')
	plt.gcf().clear()

def create_database():
	database = sqlite3.connect("gas_database.db")
	cursor = database.cursor()

	tablas = [
			"""
				CREATE TABLE IF NOT EXISTS ideal_gas(
					temperature REAL NOT NULL,
					theo_pressure REAL NOT NULL,
					exp_pressure REAL NOT NULL,
					rel_error REAL NOT NULL
				);
			"""
		]
	for tabla in tablas:
		cursor.execute(tabla);

	tablas = [
			"""
				CREATE TABLE IF NOT EXISTS real_gas(
					temperature REAL NOT NULL,
					theo_pressure REAL NOT NULL,
					exp_pressure REAL NOT NULL,
					rel_error_pressure REAL NOT NULL,
					theo_avg_vel REAL NOT NULL,
					exp_avg_vel REAL NOT NULL,
					rel_error_vel REAL NOT NULL
				);
			"""
		]
	for tabla in tablas:
		cursor.execute(tabla);

	database.commit()

def fill_database(name, table, interrogants, values):
	database = sqlite3.connect('%s.db' % name)
	cursor = database.cursor()

	sentencia = "INSERT INTO %s VALUES %s" % (table, interrogants) 
	cursor.execute(sentencia, values)
	database.commit()
	
def read_database(name, table, name_2, order):
	database = sqlite3.connect('%s.db' % name)
	cursor = database.cursor()
	sentencia = "SELECT * FROM '%s' ORDER BY %s %s;" % (table, name_2, order)
 
	cursor.execute(sentencia)
	
	lines = cursor.fetchall()

	cursor.close()
	
	return lines

###########################################################
#														  #
#						CLASSES						  	  #
#													  	  #
###########################################################

class IdealGas(object):
	def __init__(self, N, L, T, m, R, dt, T_MAX):
		self.N = N #Number of particles
		self.L = L #Length of the box
		self.T = T #Temperature of the system
		self.m = m #Particle mass (Oxygen  2.66e-23) (He 6.64e-24)
		self.R = R #Particle radius
		self.dt = dt #Step size
		self.T_MAX = T_MAX #Maximum simulation time

		self.Nt = int(self.T_MAX / self.dt) #Number of simulation steps/iterations

		self.r = np.zeros((self.N, 3)) #Array of vector positions
		self.v = np.zeros((self.N, 3)) #Array of vector velocities
		self.p = np.zeros(int(self.Nt)) # Array of momentum exchanged with the wall
		self.CM = (0,0,0) #Center of masses

		self.CM_t = []
		
	def init_particles(self):

		self.r = np.random.rand(self.N, 3) * 2 * (self.L/2 - self.R) - (self.L/2 - self.R)

		#vmed = np.sqrt(8 * k_B * T / (np.pi * m))
		vrms = np.sqrt(3*k_B * self.T / self.m)

		v_polar = np.random.random((self.N, 2))

		self.v[:,0] = (np.sin(v_polar[:,0] * np.pi) * np.cos(v_polar[:,1] * 2 * np.pi)) 
		self.v[:,1] = (np.sin(v_polar[:,0] * np.pi) * np.sin(v_polar[:,1] * 2 * np.pi)) 
		self.v[:,2] = np.cos(v_polar[:,0] * np.pi)

		self.v *= vrms

	def check_walls(self, i, x): #Checks for particle-wall collisions and return the exchanged momentum
		pex = 0
		for k in range(3):
			if abs(self.r[i][k]) + self.R > x/2:
				self.v[i][k] *= -1
				if self.r[i][k] < 0:
					self.r[i][k] = -x/2 + self.R
				else :
					self.r[i][k] = x/2 - self.R

				pex += 2 * self.m * abs(self.v[i][k])
		return pex

	def simulation(self):
		self.init_particles()

		for k in range(self.Nt):
			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]

			self.CM_t.append(np.sum(self.r, axis=0) / self.N)

		return self.p

	def isothermal_expansion(self, pull_period, L0, dL):
		self.init_particles()

		Vs = []
		Vs_discrete = []
		theo_pr_discrete = []
		time = []
		Exp_pr = []
		Theo_pr = []


		Vs_discrete.append(L0**3)	
		theo_pr_discrete.append(self.N * k_B * self.T / (L0**3))

		self.L = L0

		for k in range(self.Nt):
			if k % (pull_period/self.dt) == 0 and k > 0:
				self.L = self.L + dL
				Vs_discrete.append(self.L**3)
				theo_pr_discrete.append(self.N * k_B* self.T/(self.L**3))

			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]

			if k > 1:

				Exp_P = self.p[k] / (self.dt * 6 * self.L**2)
				Theo_P = self.N * k_B * self.T / (self.L**3)

				Exp_pr.append(Exp_P)
				Theo_pr.append(Theo_P)
				Vs.append(self.L**3)
				time.append(k*self.dt)

		if pull_period > self.dt :		
			a = int(self.T_MAX/pull_period) #Times the pressure will remain constant (as V is doesn't change)

			average_t = [] #average of the pressure in this times

			for h in range(a):
				avrg = 0
				realtime = int(pull_period/self.dt)
				for s in range(h*realtime, (h+1)*realtime-2):
					avrg = avrg + Exp_pr[s]
				average_t.append(avrg/realtime)

		plt.plot(Vs, Theo_pr, ':o', label = 'Theoretical pressure')
		plt.plot(Vs, Exp_pr, 'x', label = 'Experimental pressure')

		if pull_period > self.dt :
			plt.plot(Vs_discrete, average_t, ls = 'none', marker = 'o', label = 'Average experimental pressure')

		plt.title('Isothermal compression (T=%.2f K)' % self.T)
		plt.xlabel('Volume ($m^3$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()
		
		os.chdir(carpeta_data)
		plt.savefig('Expansion_volume.png')
		plt.gcf().clear()

		plt.plot(time, Theo_pr, '-', label = 'Theoretical pressure')
		plt.plot(time, Exp_pr, 'x', label = 'Experimental pressure')

		if pull_period > self.dt :
			plt.plot(np.linspace(0+ pull_period/2, self.T_MAX - pull_period/2, (a)), average_t, ls = 'none', marker='o', label = 'Average experimental pressure')

		plt.title('Isothermal compression (T=%.2f K)' % self.T)
		plt.xlabel('time ($s$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()
		
		os.chdir(carpeta_data)
		plt.savefig('Expansion_time.png')
		plt.gcf().clear()			

	def isochoric_heating(self, heat_period, T0, dT):	
		T_discrete = []
		T_discrete.append(T0)
		counter =0

		Exp_pr = []
		Theo_pr = []
		temperatures = []
		time = []

		self.init_particles()

		for k in range(self.Nt):

			if k % (heat_period / self.dt) == 0 and k > 0:
				self.T = self.T + dT
				T_discrete.append(self.T)

				vmed = np.sqrt(8 * k_B * self.T / (np.pi * self.m)) #Re-init particles velocities for each T
				v_polar = np.random.random((self.N, 2))

				self.v[:,0] = (np.sin(v_polar[:,0] * np.pi) * np.cos(v_polar[:,1] * 2 * np.pi)) 
				self.v[:,1] = (np.sin(v_polar[:,0] * np.pi) * np.sin(v_polar[:,1] * 2 * np.pi)) 
				self.v[:,2] = np.cos(v_polar[:,0] * np.pi)

				self.v *= vmed

			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]

			Exp_P = self.p[k]/(self.dt * 6 * self.L**2)
			Theo_P = self.N * k_B * self.T / (self.L**3)

			Exp_pr.append(Exp_P)
			Theo_pr.append(Theo_P)
			temperatures.append(self.T)
			time.append(k * self.dt)

		if heat_period > self.dt:
			a = int(self.T_MAX / heat_period) #Times the pressure will remain constant (as T doesn't change)

			average_t = [] #average of the pressure in this times

			for h in range(a):
				avrg = 0
				realtime = int(heat_period / self.dt)
				for s in range(h*realtime, (h+1)*realtime-2):
					avrg = avrg + Exp_pr[s]
				average_t.append(avrg/realtime)

		plt.plot(temperatures, Theo_pr, ':.', label = 'Theoretical pressure')
		plt.plot(temperatures, Exp_pr, 'x', label = 'Experimental pressure')

		if heat_period > self.dt:
			plt.plot(T_discrete, average_t, '-o', label = 'Average experimental pressure')

		plt.title('Isochoric process (V=%.2f $m^3$)'%(self.L**3))
		plt.xlabel('Temperature ($K$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()

		os.chdir(carpeta_data)
		plt.savefig('Heating_temperature.png')
		plt.gcf().clear()

		plt.plot(time, Theo_pr, '-', label = 'Theoretical pressure')
		plt.plot(time, Exp_pr, 'x', label = 'Experimental pressure')

		if heat_period > self.dt:
			plt.plot(np.linspace(0+ heat_period/2, self.T_MAX - heat_period/2, (a)), average_t, ls = 'none', marker='o', label = 'Average experimental pressure')

		plt.title('Isochoric process (V=%.2f $m^3$)'%(self.L**3))
		plt.xlabel('time ($s$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()

		os.chdir(carpeta_data)
		plt.savefig('Heating_time.png')
		plt.gcf().clear()

	def animation(self):

		self.init_particles()

		def update(t, lines):
			k = int(t / self.dt)

			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]

			self.CM = np.sum(self.r, axis=0) / self.N
			self.CM_t.append(self.CM)

			lines[0].set_data(self.r[:,0], self.r[:,1])
			lines[0].set_3d_properties(self.r[:,2])

			lines[1].set_data([self.CM[0]], [self.CM[1]])
			lines[1].set_3d_properties([self.CM[2]])

			return lines

		# Attaching 3D axis to the figure
		fig = plt.figure()
		ax = p3.Axes3D(fig)


		# Setting the axes properties
		ax.set_xlim3d([-self.L/2, self.L/2])
		ax.set_xlabel('X')

		ax.set_ylim3d([-self.L/2, self.L/2])
		ax.set_ylabel('Y')

		ax.set_zlim3d([-self.L/2, self.L/2])
		ax.set_zlabel('Z')

		T1 = calculate_T(self.v, self.m, self.N)

		lines = []
		lines.append(ax.plot(self.r[:,0], self.r[:,1], self.r[:,2], ls='None', marker='.', label = 'Particles at $T_0=%.2f\,K$'%T1)[0])
		lines.append(ax.plot([self.CM[0]], [self.CM[1]], [self.CM[2]], marker='o', color='r', label='Center of masses')[0])

		ani = FuncAnimation(fig, update, fargs=(lines,), frames=np.linspace(0, self.T_MAX-self.dt, int(self.T_MAX / self.dt)),
		                    blit=True, interval=1, repeat=False)
		#ani.save('animation.mp4', fps=20, writer="ffmpeg", codec="libx264")
		plt.legend(loc = 'upper left')
		plt.show()

	def experiments(self):

		Exp_P, Theo_P = calculate_P(self.N, self.L, self.T, self.dt, self.p)
		T1 = calculate_T(self.v, self.m, self.N)

		A = 6*self.L**2

		theo_P_2 = self.N * k_B * T1 / (self.L**3)

		rel_error = abs(Theo_P - Exp_P) / Theo_P * 100
		rel_error2 = abs(theo_P_2 - Exp_P) / theo_P_2 * 100

		os.chdir(carpeta_data)
		fill_database('gas_database', 'ideal_gas', '(?,?,?,?)', [T1, Theo_P, Exp_P, rel_error])

		#Pressure plot
		theo_array = []
		exp_array = []

		for i in range(self.Nt):
			theo_array.append(Theo_P)
			exp_array.append(Exp_P)

		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), self.p / (self.dt * A), ':', label = 'Experimental pressure')
		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), theo_array, label = 'Theoretical pressure')
		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), exp_array, label = 'Average experimental pressure')
		plt.title('Pressure of the system')
		plt.xlabel('t(s)')
		plt.ylabel('P(Pa)')
		plt.legend()

		os.chdir(carpeta_data)
		plt.savefig('Pressure.png')
		plt.gcf().clear()

		#CM PLOTS
		transposed = np.transpose(self.CM_t)
		t = np.linspace(0, self.T_MAX - self.dt, len(self.CM_t))

		std_x = np.std(transposed[:][0])
		std_y = np.std(transposed[:][1])
		std_z = np.std(transposed[:][2])

		mean_x = np.mean(transposed[:][0])
		mean_y = np.mean(transposed[:][1])
		mean_z = np.mean(transposed[:][2])

		plt.subplot(2, 2, 1)
		plt.plot(t, transposed[:][0]) #CM_x
		meanx = plt.plot(t, np.linspace(mean_x, mean_x, len(self.CM_t)), color='red', label='Mean') #mean
		sdx = plt.plot(t, np.linspace(std_x, std_x, len(self.CM_t)), ls = ':', color='green', label=r'$\pm$ SD') #+ standard deviation
		plt.plot(t, np.linspace(-std_x, -std_x, len(self.CM_t)), ls = ':', color='green') #- standard deviation
		plt.xlabel('t')
		plt.ylabel('x')

		plt.subplot(2, 2, 2)
		plt.plot(t, transposed[:][1]) #CM_y
		plt.plot(t, np.linspace(mean_y, mean_y, len(self.CM_t)), color='red', label='Mean') #mean
		plt.plot(t, np.linspace(std_y, std_y, len(self.CM_t)), ls = ':', color='green', label=r'$\pm$ SD') #+ standard deviation
		plt.plot(t, np.linspace(-std_y, -std_y, len(self.CM_t)), ls = ':', color='green') #- standard deviation
		plt.xlabel('t')
		plt.ylabel('y')

		plt.subplot(2, 2, 3)
		plt.plot(t, transposed[:][2]) #CM_z
		plt.plot(t, np.linspace(mean_z, mean_z, len(self.CM_t)), color='red', label='Mean') #mean
		plt.plot(t, np.linspace(std_z, std_z, len(self.CM_t)), ls = ':', color='green', label=r'$\pm$ SD') #+ standard deviation
		plt.plot(t, np.linspace(-std_z, -std_z, len(self.CM_t)), ls = ':', color='green') #- standard deviation
		plt.xlabel('t')
		plt.ylabel('z')

		plt.suptitle(r'$\vec{r}_{CM}$ vs t')
		plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
		
		os.chdir(carpeta_data)
		plt.savefig('CM_plots.png')
		plt.gcf().clear()

class RealGas(object):
	def __init__(self, N, L, T, m, R, dt, T_MAX):
		self.N = N #Number of particles
		self.L = L #Length of the box
		self.T = T #Temperature of the system
		self.m = m #Particle mass (Oxygen  2.66e-23) (He 6.64e-24)
		self.R = R #Particle radius
		self.dt = dt #Step size
		self.T_MAX = T_MAX #Maximum simulation time

		self.Nt = int(self.T_MAX / self.dt) #Number of simulation steps/iterations

		self.r = np.zeros((self.N, 3)) #Array of vector positions
		self.v = np.zeros((self.N, 3)) #Array of vector velocities
		self.p = np.zeros(int(self.Nt)) # Array of momentum exchanged with the wall
		self.CM = (0,0,0) #Center of masses

		self.hist = []
		self.CM_t = []

	def init_particles(self):

		self.r = np.random.rand(self.N, 3) * 2 * (self.L/2 - self.R) - (self.L/2 - self.R)

		#vmed = np.sqrt(8 * k_B * T / (np.pi * m))
		vrms = np.sqrt(3*k_B * self.T / self.m)

		v_polar = np.random.random((self.N, 2))

		self.v[:,0] = (np.sin(v_polar[:,0] * np.pi) * np.cos(v_polar[:,1] * 2 * np.pi)) 
		self.v[:,1] = (np.sin(v_polar[:,0] * np.pi) * np.sin(v_polar[:,1] * 2 * np.pi)) 
		self.v[:,2] = np.cos(v_polar[:,0] * np.pi)

		self.v *= vrms

	def check_walls(self, i, x): #Checks for particle-wall collisions and return the exchanged momentum
		pex = 0
		for k in range(3):
			if abs(self.r[i][k]) + self.R > x/2:
				self.v[i][k] *= -1
				if self.r[i][k] < 0:
					self.r[i][k] = -x/2 + self.R
				else :
					self.r[i][k] = x/2 - self.R

				pex += 2 * self.m * abs(self.v[i][k])
		return pex

	def check_collision(self):  #Check for collisions with other particles of te same gas

		dists = np.sqrt(mod(self.r - self.r[:,np.newaxis])) 
		cols2 = (0 < dists) & (dists < 2*self.R)
		idx_i, idx_j = np.nonzero(cols2)
		# ***possibility to simplify this *** #
		for i, j in zip(idx_i, idx_j):
		    if j < i:
		        # skip duplications and same particle
		        continue 

		    rij = self.r[i] - self.r[j]
		    d = mod(rij)
		    vij = self.v[i] - self.v[j]
		    dv = np.dot(vij, rij) * rij / d
		    self.v[i] -= dv
		    self.v[j] += dv
			
		    # update the positions so they are no longer in contact
		    self.r[i] += self.dt * self.v[i]
		    self.r[j] += self.dt * self.v[j]

	def simulation(self):
		self.init_particles()

		for k in range(self.Nt):
			self.check_collision()

			modv = 0
			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]
				modv = modv + LA.norm(self.v[i])

			self.CM_t.append(np.sum(self.r, axis=0) / self.N)

			average_v = modv / self.N
			self.hist.append(average_v)

		return self.p

	def isothermal_expansion(self, pull_period, L0, dL):
		self.init_particles()

		Vs = []
		Vs_discrete = []
		theo_pr_discrete = []
		time = []
		Exp_pr = []
		Theo_pr = []


		Vs_discrete.append(L0**3)	
		theo_pr_discrete.append(self.N * k_B * self.T / (L0**3))

		self.L = L0

		for k in range(self.Nt):
			self.check_collision()

			if k % (pull_period/self.dt) == 0 and k > 0:
				self.L = self.L + dL
				Vs_discrete.append(self.L**3)
				theo_pr_discrete.append(self.N * k_B* self.T/(self.L**3))

			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]

			if k > 1:

				Exp_P = self.p[k] / (self.dt * 6 * self.L**2)
				Theo_P = self.N * k_B * self.T / (self.L**3)

				Exp_pr.append(Exp_P)
				Theo_pr.append(Theo_P)
				Vs.append(self.L**3)
				time.append(k*self.dt)

		if pull_period > self.dt :		
			a = int(self.T_MAX/pull_period) #Times the pressure will remain constant (as V is doesn't change)

			average_t = [] #average of the pressure in this times

			for h in range(a):
				avrg = 0
				realtime = int(pull_period/self.dt)
				for s in range(h*realtime, (h+1)*realtime-2):
					avrg = avrg + Exp_pr[s]
				average_t.append(avrg/realtime)

		plt.plot(Vs, Theo_pr, ':o', label = 'Theoretical pressure')
		plt.plot(Vs, Exp_pr, 'x', label = 'Experimental pressure')

		if pull_period > self.dt :
			plt.plot(Vs_discrete, average_t, ls = 'none', marker = 'o', label = 'Average experimental pressure')

		plt.title('Isothermal compression (T=%.2f K)' % self.T)
		plt.xlabel('Volume ($m^3$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()
		
		os.chdir(carpeta_data)
		plt.savefig('Expansion_volume.png')
		plt.gcf().clear()

		plt.plot(time, Theo_pr, '-', label = 'Theoretical pressure')
		plt.plot(time, Exp_pr, 'x', label = 'Experimental pressure')

		if pull_period > self.dt :
			plt.plot(np.linspace(0+ pull_period/2, self.T_MAX - pull_period/2, (a)), average_t, ls = 'none', marker='o', label = 'Average experimental pressure')

		plt.title('Isothermal compression (T=%.2f K)' % self.T)
		plt.xlabel('time ($s$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()
		
		os.chdir(carpeta_data)
		plt.savefig('Expansion_time.png')
		plt.gcf().clear()			

	def isochoric_heating(self, heat_period, T0, dT):	
		T_discrete = []
		T_discrete.append(T0)
		counter =0

		Exp_pr = []
		Theo_pr = []
		temperatures = []
		time = []

		self.init_particles()

		for k in range(self.Nt):
			self.check_collision()

			if k % (heat_period / self.dt) == 0 and k > 0:
				self.T = self.T + dT
				T_discrete.append(self.T)

				vmed = np.sqrt(8 * k_B * self.T / (np.pi * self.m)) #Re-init particles velocities for each T
				v_polar = np.random.random((self.N, 2))

				self.v[:,0] = (np.sin(v_polar[:,0] * np.pi) * np.cos(v_polar[:,1] * 2 * np.pi)) 
				self.v[:,1] = (np.sin(v_polar[:,0] * np.pi) * np.sin(v_polar[:,1] * 2 * np.pi)) 
				self.v[:,2] = np.cos(v_polar[:,0] * np.pi)

				self.v *= vmed

			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]

			Exp_P = self.p[k]/(self.dt * 6 * self.L**2)
			Theo_P = self.N * k_B * self.T / (self.L**3)

			Exp_pr.append(Exp_P)
			Theo_pr.append(Theo_P)
			temperatures.append(self.T)
			time.append(k * self.dt)

		if heat_period > self.dt:
			a = int(self.T_MAX / heat_period) #Times the pressure will remain constant (as T doesn't change)

			average_t = [] #average of the pressure in this times

			for h in range(a):
				avrg = 0
				realtime = int(heat_period / self.dt)
				for s in range(h*realtime, (h+1)*realtime-2):
					avrg = avrg + Exp_pr[s]
				average_t.append(avrg/realtime)

		plt.plot(temperatures, Theo_pr, ':.', label = 'Theoretical pressure')
		plt.plot(temperatures, Exp_pr, 'x', label = 'Experimental pressure')

		if heat_period > self.dt:
			plt.plot(T_discrete, average_t, '-o', label = 'Average experimental pressure')

		plt.title('Isochoric process (V=%.2f $m^3$)'%(self.L**3))
		plt.xlabel('Temperature ($K$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()

		os.chdir(carpeta_data)
		plt.savefig('Heating_temperature.png')
		plt.gcf().clear()

		plt.plot(time, Theo_pr, '-', label = 'Theoretical pressure')
		plt.plot(time, Exp_pr, 'x', label = 'Experimental pressure')

		if heat_period > self.dt:
			plt.plot(np.linspace(0+ heat_period/2, self.T_MAX - heat_period/2, (a)), average_t, ls = 'none', marker='o', label = 'Average experimental pressure')

		plt.title('Isochoric process (V=%.2f $m^3$)'%(self.L**3))
		plt.xlabel('time ($s$)')
		plt.ylabel('Pressure (Pa)')
		plt.legend()

		os.chdir(carpeta_data)
		plt.savefig('Heating_time.png')
		plt.gcf().clear()


	def animation(self):

		self.init_particles()

		def update(t, lines):
			k = int(t / self.dt)
			self.check_collision(self.r, self.v)

			modv = 0
			for i in range(self.N):
				self.p[k] += self.check_walls(i, self.L)
				self.r[i] += self.dt * self.v[i]
				modv = modv + LA.norm(self.v[i])

			self.CM = np.sum(self.r, axis=0) / self.N
			self.CM_t.append(self.CM)

			average_v = modv / self.N
			self.hist.append(average_v)

			lines[0].set_data(self.r[:,0], self.r[:,1])
			lines[0].set_3d_properties(self.r[:,2])

			lines[1].set_data([self.CM[0]], [self.CM[1]])
			lines[1].set_3d_properties([self.CM[2]])

			return lines

		# Attaching 3D axis to the figure
		fig = plt.figure()
		ax = p3.Axes3D(fig)


		# Setting the axes properties
		ax.set_xlim3d([-self.L/2, self.L/2])
		ax.set_xlabel('X')

		ax.set_ylim3d([-self.L/2, self.L/2])
		ax.set_ylabel('Y')

		ax.set_zlim3d([-self.L/2, self.L/2])
		ax.set_zlabel('Z')

		T1 = calculate_T(self.v, self.m, self.N)

		lines = []
		lines.append(ax.plot(self.r[:,0], self.r[:,1], self.r[:,2], ls='None', marker='.', label = 'Particles at $T_0=%.2f\,K$'%T1)[0])
		lines.append(ax.plot([self.CM[0]], [self.CM[1]], [self.CM[2]], marker='o', color='r', label='Center of masses')[0])

		ani = FuncAnimation(fig, update, fargs=(lines,), frames=np.linspace(0, self.T_MAX-self.dt, self.Nt),
		                    blit=True, interval=1, repeat=False)
		#ani.save('animation.mp4', fps=20, writer="ffmpeg", codec="libx264")
		plt.legend(loc = 'upper left')
		plt.show()

	def experiments(self):

		Exp_P, Theo_P = calculate_P(self.N, self.L, self.T, self.dt, self.p)
		T1 = calculate_T(self.v, self.m, self.N)

		A = 6*self.L**2

		theo_P_2 = self.N * k_B * T1 / (self.L**3)

		rel_error = abs(Theo_P - Exp_P) / Theo_P * 100
		rel_error2 = abs(theo_P_2 - Exp_P) / theo_P_2 * 100

		#PRESSURE PLOT
		theo_array = []
		exp_array = []

		for i in range(self.Nt):
			theo_array.append(Theo_P)
			exp_array.append(Exp_P)

		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), self.p / (self.dt * A), ':', label = 'Experimental pressure')
		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), theo_array, label = 'Theoretical pressure')
		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), exp_array, label = 'Average experimental pressure')
		plt.title('Pressure of the system')
		plt.xlabel('t(s)')
		plt.ylabel('P(Pa)')
		plt.legend()
		
		os.chdir(carpeta_data)
		plt.savefig('Pressure.png')
		plt.gcf().clear()

		#VELOCITY PLOT
		theo_array = []
		theo_array_rms = []

		vmed = np.sqrt(8 * k_B * self.T / (np.pi * self.m))
		vrms = np.sqrt(3 * k_B * self.T / self.m)

		exp_avg_vel = np.sum(self.hist) / len(self.hist)

		rel_error_vel = abs(vmed-exp_avg_vel)/vmed * 100

		fill_database('gas_database', 'real_gas', '(?,?,?,?,?,?,?)', [T1, Theo_P, Exp_P, rel_error, vmed, exp_avg_vel, rel_error_vel])

		plt.plot(np.linspace(0, self.T_MAX - self.dt, len(self.hist)), self.hist, '-', label = 'Experimental data')
		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), np.linspace(vmed, vmed, self.Nt), label = r'Maxwell-Boltzman $v_{med}$')
		plt.plot(np.linspace(0, self.T_MAX - self.dt, self.Nt), np.linspace(vrms, vrms, self.Nt), label = r'Maxwell-Boltzman $v_{rms}$')
		plt.title('Average velocity ($T=%.2f\,K$)' % T1)
		plt.xlabel('t(s)')
		plt.ylabel('$v_{med}$(m/s)')
		plt.legend()
		
		os.chdir(carpeta_data)
		plt.savefig('Velocity.png')
		plt.gcf().clear()

		#VELOCITY HISTOGRAM

		modv = []
		for i in range(len(self.v)):
			mod = self.v[i][0]**2 + self.v[i][1]**2 + self.v[i][2]**2
			modv.append(np.sqrt(mod))

		v_max = 3*np.sqrt((3*k_B* self.T)/ self.m)
		dv = 0.01

		vel_hist(modv, self.v, v_max, dv, self.m, self.T)

		#CENTER OF MASSES MOVEMENT PLOTS

		transposed = np.transpose(self.CM_t)
		t = np.linspace(0, self.T_MAX - self.dt, len(self.CM_t))

		std_x = np.std(transposed[:][0])
		std_y = np.std(transposed[:][1])
		std_z = np.std(transposed[:][2])

		mean_x = np.mean(transposed[:][0])
		mean_y = np.mean(transposed[:][1])
		mean_z = np.mean(transposed[:][2])

		plt.subplot(2, 2, 1)
		plt.plot(t, transposed[:][0]) #CM_x
		meanx = plt.plot(t, np.linspace(mean_x, mean_x, len(self.CM_t)), color='red', label='Mean') #mean
		sdx = plt.plot(t, np.linspace(std_x, std_x, len(self.CM_t)), ls = ':', color='green', label=r'$\pm$ SD') #+ standard deviation
		plt.plot(t, np.linspace(-std_x, -std_x, len(self.CM_t)), ls = ':', color='green') #- standard deviation
		plt.xlabel('t')
		plt.ylabel('x')

		plt.subplot(2, 2, 2)
		plt.plot(t, transposed[:][1]) #CM_y
		plt.plot(t, np.linspace(mean_y, mean_y, len(self.CM_t)), color='red', label='Mean') #mean
		plt.plot(t, np.linspace(std_y, std_y, len(self.CM_t)), ls = ':', color='green', label=r'$\pm$ SD') #+ standard deviation
		plt.plot(t, np.linspace(-std_y, -std_y, len(self.CM_t)), ls = ':', color='green') #- standard deviation
		plt.xlabel('t')
		plt.ylabel('y')

		plt.subplot(2, 2, 3)
		plt.plot(t, transposed[:][2]) #CM_z
		plt.plot(t, np.linspace(mean_z, mean_z, len(self.CM_t)), color='red', label='Mean') #mean
		plt.plot(t, np.linspace(std_z, std_z, len(self.CM_t)), ls = ':', color='green', label=r'$\pm$ SD') #+ standard deviation
		plt.plot(t, np.linspace(-std_z, -std_z, len(self.CM_t)), ls = ':', color='green') #- standard deviation
		plt.xlabel('t')
		plt.ylabel('z')

		plt.suptitle(r'$\vec{r}_{CM}$ vs t')
		plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')

		os.chdir(carpeta_data)
		plt.savefig('CM_plots.png')
		plt.gcf().clear()


#Fer isothermal and isochoric processes a la classe de la interfase per utilitzar la barra de carrega
#Boltzmann constant (J/K) ########################################################### #														  # #						FUNCTIONS						  # #													  	  # ########################################################### Computes the squared sum over the last axis of the numpy.ndarray v #Maxwell-Boltzman distribution CREATE TABLE IF NOT EXISTS ideal_gas(
 					temperature REAL NOT NULL,
 					theo_pressure REAL NOT NULL,
 					exp_pressure REAL NOT NULL,
 					rel_error REAL NOT NULL
 				); CREATE TABLE IF NOT EXISTS real_gas(
 					temperature REAL NOT NULL,
 					theo_pressure REAL NOT NULL,
 					exp_pressure REAL NOT NULL,
 					rel_error_pressure REAL NOT NULL,
 					theo_avg_vel REAL NOT NULL,
 					exp_avg_vel REAL NOT NULL,
 					rel_error_vel REAL NOT NULL
 				); ########################################################### #														  # #						CLASSES						  	  # #													  	  # ########################################################### #Number of particles #Length of the box #Temperature of the system #Particle mass (Oxygen  2.66e-23) (He 6.64e-24) #Particle radius #Step size #Maximum simulation time #Number of simulation steps/iterations #Array of vector positions #Array of vector velocities # Array of momentum exchanged with the wall #Center of masses #vmed = np.sqrt(8 * k_B * T / (np.pi * m)) #Checks for particle-wall collisions and return the exchanged momentum #Times the pressure will remain constant (as V is doesn't change) #average of the pressure in this times #Re-init particles velocities for each T #Times the pressure will remain constant (as T doesn't change) #average of the pressure in this times # Attaching 3D axis to the figure # Setting the axes properties #ani.save('animation.mp4', fps=20, writer="ffmpeg", codec="libx264") #Pressure plot #CM PLOTS #CM_x #mean #+ standard deviation #- standard deviation #CM_y #mean #+ standard deviation #- standard deviation #CM_z #mean #+ standard deviation #- standard deviation #Number of particles #Length of the box #Temperature of the system #Particle mass (Oxygen  2.66e-23) (He 6.64e-24) #Particle radius #Step size #Maximum simulation time #Number of simulation steps/iterations #Array of vector positions #Array of vector velocities # Array of momentum exchanged with the wall #Center of masses #vmed = np.sqrt(8 * k_B * T / (np.pi * m)) #Checks for particle-wall collisions and return the exchanged momentum #Check for collisions with other particles of te same gas # ***possibility to simplify this *** # # skip duplications and same particle # update the positions so they are no longer in contact #Times the pressure will remain constant (as V is doesn't change) #average of the pressure in this times #Re-init particles velocities for each T #Times the pressure will remain constant (as T doesn't change) #average of the pressure in this times # Attaching 3D axis to the figure # Setting the axes properties #ani.save('animation.mp4', fps=20, writer="ffmpeg", codec="libx264") #PRESSURE PLOT #VELOCITY PLOT #VELOCITY HISTOGRAM #CENTER OF MASSES MOVEMENT PLOTS #CM_x #mean #+ standard deviation #- standard deviation #CM_y #mean #+ standard deviation #- standard deviation #CM_z #mean #+ standard deviation #- standard deviation #Fer isothermal and isochoric processes a la classe de la interfase per utilitzar la barra de carrega