...
 
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -7,41 +7,44 @@ def chi(A,x):
function to determine chi_A(t)
last modified: 20/03/19
input:
A: input data
(1,n_max)-array
x: x data
(1,n_max-1)-array
output:
chi: autocorrelation function
(1,x)-array
Parameters
----------
A : (1,n_max)-array
input data
x : (1,n_max-1)-array
x data
Returns
-------
chi : (1,x)-array
autocorrelation function
"""
# Process input
n_max = np.size(A)
A_tile = np.tile(A,[n_max,1]) # (n_max,n_max)-array
A_tile = np.tile(A,[n_max,1]) # (n_max,n_max)-array
# Calculate scalefactor N
N = 1/(n_max-x) # (1,x)-array
N = 1/(n_max-x) # (1,x)-array
# Calculate chi
# # chi_1(t) = 1/(n_max-t) * sum_{n=0}^{n_max-t} A(n)A(n+t)
AA = np.convolve(A,np.flip(A,0))[0:n_max][::-1] # sum_{n=0}^{n_max-t} A(n)A(n+t)
chi_1 = np.multiply(N,AA) # (1,n_max)-array
AA = np.convolve(A,np.flip(A,0))[0:n_max][::-1]
chi_1 = np.multiply(N,AA) # (1,n_max)-array
# # chi_2(t) = 1/(n_max-t) * sum_{n=0}^{n_max-t} A(n)
A_triu_2 = np.triu(np.flip(A_tile,1))
A_2 = np.sum(A_triu_2,1)
chi_2 = np.multiply(N,A_2) # (1,n_max)-array
chi_2 = np.multiply(N,A_2) # (1,n_max)-array
# # chi_3(t) = 1/(n_max-t) * sum_{n=0}^{n_max-t} A(n+1)
A_triu_3 = np.triu(A_tile)
A_3 = np.sum(A_triu_3,1)
chi_3 = np.multiply(N,A_3) # (1,n_max)-array
chi_3 = np.multiply(N,A_3) # (1,n_max)-array
# # chi = chi_1 + chi_2*chi_3
chi = chi_1 - np.multiply(chi_2,chi_3) # (1,n_max)-array
chi = chi_1 - np.multiply(chi_2,chi_3) # (1,n_max)-array
# return
return chi
......@@ -51,17 +54,21 @@ def test_func(x, a, b):
exponential function to fit chi
last modified: 26/03/2019
input:
x: x_data
(1,N) array
a: amplitude
(float)
b: tau
(float)
output:
a*np.exp(-x/b)+c: function
(1,N) array
Parameters
----------
x : (1,N) array
x_data
a : (float)
amplitude
b : (float)
tau
Returns
-------
a*np.exp(-x/b) : (1,N) array
function
"""
return a*np.exp(-x/b)
......@@ -71,13 +78,16 @@ def tau(y_data):
function to find correlation time of the simulation (tau) by fitting chi to exponential test function
last modified: 26/03/2019
input:
y_data: data to which function is fitted
(1,N) array
Parameters
----------
y_data : (1,N) array
data to which function is fitted
output:
tau: correlation time of the simulation
(float)
Returns
-------
tau : (float)
correlation time of the simulation
"""
# data
x_data = np.array(range(len(y_data)))
......@@ -93,17 +103,20 @@ def calc_std_mean(tau,data):
function to calculate standard deviation (std) and mean (mean) of data using data blocking
last modified: 20/03/19
input:
tau: correlation time of the simulation
(float)
data: data from which std and mean are calculated
(1,N)-array
Parameters
----------
tau : float
correlation time of the simulation
data : (1,N)-array
data from which std and mean are calculated
output:
data_std: standard deviation of data
(float)
data_mean: mean of data
(float)
Returns
-------
data_std : float
standard deviation of data
data_mean : float
mean of data
"""
# define length of data block
......@@ -131,16 +144,17 @@ def error_combine(data):
overhead function to calculate standard deviation and mean using datablocking
last modified: 20/03/2019
input:
data: data from which std and mean are calculated
(1,N)-array
output:
data_std: standard deviation of data
(float)
data_mean: mean of data
(float)
Parameters
----------
data : (1,N)-array
data from which std and mean are calculated
Returns
-------
data_std : (float)
standard deviation of data
data_mean : (float)
mean of data
"""
# process input
n_max = np.size(data)
......
......@@ -19,19 +19,32 @@ def energy(K, U, std_E, std_K, std_U, error_alpha, h, N, units='physical'):
"""
Plots kinetic, potential and total energy as function of time.
Input:
K Total kinetic energy of the system for every step ([1,N_steps - 1] array)
U Total potential energy of the system for every step ([1,N_steps - 1] array)
std_E Standard deviation of the total energy
std_K Standard deviation of the kinetic energy
std_U Standard deviation of the potential energy
error_alpha Transparancy of errorbar shade
h Length of time step
N Number of particles
units Natural or physical units as output
Parameters
----------
K : ([1,N_steps - 1] array)
Total kinetic energy of the system for every step
U : ([1,N_steps - 1] array)
Total potential energy of the system for every step
std_E : float
Standard deviation of the total energy
std_K : float
Standard deviation of the kinetic energy
std_U : float
Standard deviation of the potential energy
error_alpha : float
Transparancy of errorbar shade
h : float
Length of time step
N : float
Number of particles
units : string
Natural or physical units as output
Output:
Raises
------
plot
IOError
When units field unequel to either 'natural' or 'physical'
"""
if units=='natural':
......@@ -78,18 +91,30 @@ def pressure(P, std_P, error_alpha, h, density, K, N, units='physical'):
"""
Plots the pressure inside the particle box as a function of time.
Input:
P Array with the pressure inside the particle box at each time step ((N_time_steps) np.array)
std_P Standard deviation of the pressure
error_alpha Transparancy of errorbar shade
h Length of time step
density Particle density inside the box
K Total kinetic energy of the system for every step ([1,N_steps - 1] array)
N Number of particles
units Natural or physical units as output
Parameters
----------
P : ((N_time_steps) np.array)
Array with the pressure inside the particle box at each time step
std_P : float
Standard deviation of the pressure
error_alpha : float
Transparancy of errorbar shade
h : float
Length of time step
density : float
Particle density inside the box
K : ([1,N_steps - 1] array)
Total kinetic energy of the system for every step
N : float
Number of particles
units : string
Natural or physical units as output
Output:
Raises
------
plot
IOError
When units field unequel to either 'natural' or 'physical'
"""
if units=='natural':
......@@ -126,17 +151,28 @@ def specific_heat(CV, std_CV, error_alpha, N_CV, h, N, units='physical'):
"""
Plots the specific heat (heat capacity) as a function of time.
Input:
CV Specific heat
std_CV Standard deviation of the specific heat
error_alpha Transparancy of errorbar shade
N_CV Number of timesteps to average over
h Length of time step
N Number of particles
units Natural or physical units as output
Parameters
----------
CV : array
Specific heat
std_CV : float
Standard deviation of the specific heat
error_alpha : float
Transparancy of errorbar shade
N_CV : float
Number of timesteps to average over
h : float
Length of time step
N : float
Number of particles
units : string
Natural or physical units as output
Output:
Raises
------
plot
IOError
When units field unequel to either 'natural' or 'physical'
"""
CV = CV/N
......@@ -170,14 +206,22 @@ def pair_correlation(r, G, N, units='physical'):
"""
Plots the pair correlation as a function of distance.
Input:
G Pair correlation function
r Distances for which the pair correlation function is defined
N Number of particles
units Natural or physical units as output
Parameters
----------
G : array
Pair correlation function
r : array
Distances for which the pair correlation function is defined
N : float
Number of particles
units : string
Natural or physical units as output
Output:
Raises
------
plot
IOError
When units field unequel to either 'natural' or 'physical'
"""
if units=='natural':
pass
......@@ -202,15 +246,23 @@ def diffusion(diff, N, N_steps, h, units='physical'):
"""
Plots the diffusion as a function of time.
Input:
diff Vector with the average distance traveled by the particles with respect to t=0 ((N_time_steps) np.array)
N Number of particles inside the box
N_steps Number of timesteps
h Length of timestep
units Natural or physical units as output
Parameters
----------
diff : ((N_time_steps) np.array) Vector with the average distance traveled by the particles with respect to t=0
N : float
Number of particles inside the box
N_steps : float
Number of timesteps
h : float
Length of timestep
units : string
Natural or physical units as output
Output:
Raises
------
plot
IOError
When units field unequel to either 'natural' or 'physical'
"""
if units=='natural':
......@@ -238,14 +290,21 @@ def temperature(T, h, N, units='physical'):
"""
Plots the temperature as a function of time.
Input:
T Temperature per time step
h Length of timestep
N Number of particles inside the box
units Natural or physical units as output
Parameters
T : array
Temperature per time step
h : float
Length of timestep
N : float
Number of particles inside the box
units : string
Natural or physical units as output
Output:
Raises
------
plot
IOError
When units field unequel to either 'natural' or 'physical'
"""
if units=='natural':
......
# import modules
import numpy as np
def create_table(density, T, E_kin, P, N, U, Cv):
"""
Function for a clear comparison with literature data and data as found in by the simulated model.
Parameters
---------
density : float
T : float
temperature
E_kin : array
Kinetic energy
P : float
Pressure
N : float
Number of particles
U : array
Potential energy
Cv : float
Specific heat
Returns
-------
Table with simulated data and data found in literature.
"""
print("In natural units:")
# write table
table = [
["", "Density (set)", "Temperature (set)", "Temperature", "Pressure (beta*P/rho)", "Potential energy", "Specific heat"],
["Values book (1)", 0.88, 1.0, 0.990, 2.98, -5.704, "1.5 (Ideal gas)"],
......@@ -12,12 +36,12 @@ def create_table(density, T, E_kin, P, N, U, Cv):
["Values book (3)", 0.70, 1.0, 1.014, 1.06, -4.662, "1.5 (Ideal gas)"],
["Our values", f"{density:.2e}", f"{T:.2e}", f"{np.mean(2*E_kin[-800:]/((N-1)*3)):.2e}", f"{P:.2e}", f"{np.mean(U[-800:])/N:.2e}", f"{Cv/N:.2e}"],
]
return(table)
# return
return(table)
def print_table(table):
"""Function to print table"""
longest_cols = [
(max([len(str(row[i])) for row in table]) + 3)
for i in range(len(table[0]))
......@@ -28,5 +52,6 @@ def print_table(table):
def validation_table(density, T, E_kin, P, N, U, Cv):
"""Function to set up table for comparison with literature values"""
table = create_table(density, T, E_kin, P, N, U, Cv)
print_table(table)
\ No newline at end of file