-
Notifications
You must be signed in to change notification settings - Fork 0
/
lack_reader.py
53 lines (39 loc) · 2.1 KB
/
lack_reader.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
# -*- coding: utf-8 -*-
'''Reads the outputs from my Lack model simulations so the data can be used in graphs etc'''
import numpy as np
import pandas as pd
import lack_plots as lp
import lack_functions as lf
import matplotlib.pyplot as plt
file_name = 'lack_model_output.txt'
# file_name = '..\\MAIN\\Charge and Size Distributions\\New_Output_Files\\G90_63+_Mastersizer.txt'
fit_type = 'complex' # Simple or complex depending on the shape of the output
df = pd.read_csv(file_name, sep=',', header=16) #header may need to be altered for different versions of the code's ouptut file
info_df = pd.read_csv(file_name, sep='\t', header=0, nrows=15)
# print(df.columns.values) # Column headers
# print(max(df['final high e'])) # High energy states
vel = df[['velocity(x)', 'velocity(y)', 'velocity(z)']].to_numpy()
pos = df[['position(x)', 'position(y)', 'position(z)']].to_numpy()
radii, charge, mass = df['radii'], -df['charge'], df['masses']
box_length = lf.calc_box_len(float(info_df['Parameters: '].iloc[0]), radii)
electron_surface_density = 1
diameters = radii * 2
dimentionless_charge = charge / (4 * np.pi * electron_surface_density)
lp.size_histogram(diameters, n_bins=15)
# lp.speed_histogram(vel, 10)
# lp.scatter_speeds(vel, radii)
lp.plot_3D(pos, charge, radii, box_length / 1000)
lp.ke_histogram(vel, mass, 25, dpi=200)
lp.plot_charge(diameters, dimentionless_charge, "Diameter", dpi=200)
if fit_type == 'simple':
a, b, c = [8.951186771643852e-06, 1.6, -100]
a, b, c = lf.get_fit(diameters, dimentionless_charge, initial_abc_guess=[a, b, c], minimise_op="R2")
lp.plot_fit(a, b, c, lf.round_down_to_1sf(min(diameters)), lf.round_up_to_1sf(max(diameters)))
print(f'a, b, c = {a}, {b}, {c}')
elif fit_type == 'complex':
a, b = [0.0015572432345201592, -1156.8152071504148]
a, b = lf.get_fit_complex(diameters, dimentionless_charge, initial_ad_guess=[a, b], minimise_op="R2")
lp.plot_fit_complex(a, b, lf.round_down_to_1sf(min(diameters)), lf.round_up_to_1sf(max(diameters)))
print(f'a, b = {a}, {b}')
else:
raise ValueError(f'fit_type should be "Simple" or "complex", instead: {fit_type}')