import pandas as pd
from tqdm import tqdm, notebook
import numpy as np
import scipy
import plotly.express as px
import plotly.graph_objects as go
import kaleido
from IPython.display import Image
tqdm_disabled = True # True for website, change to False for local work
sampling = True
MIN_POINTS0 = 500
DIFF_LIMIT = 0.1
# electrolyte component masses, in g
MASS_ZnCl2 = 1.42
MASS_NH4Cl = 1.07
MASS_KI = 3.33
MASS_H2O = 8.52
MASS_TriEG = 0.56
total_mass_kg = (MASS_ZnCl2 + MASS_KI + MASS_H2O + MASS_TriEG)/1000.
TOTAL_VOLUME = 11 # electrolyte volume in mL, approx, measured by taking as much electrolyte as possible up into a 12 mL syringe
MASS_TO_RESERVOIRS = 14.37 # g of electrolyte actually loaded into system, based on weighing syringe before/after loading reservoirs
# molecular weights in g/mol
density = MASS_TO_RESERVOIRS/TOTAL_VOLUME
MW_ZnCl2 = 136.315
MW_NH4Cl = 53.49
MW_KI = 166.0028
MW_H2O = 18.01528
MW_TriEG = 150.174
molality_ZnCl2 = MASS_ZnCl2/MW_ZnCl2/total_mass_kg
molality_NH4Cl = MASS_NH4Cl/MW_NH4Cl/total_mass_kg
molality_KI = MASS_ZnCl2/MW_KI/total_mass_kg
molality_TriEG = MASS_TriEG/MW_TriEG/total_mass_kg
molality_H2O = MASS_H2O/MW_H2O/total_mass_kg
molarity_ZnCl2 = MASS_ZnCl2/MW_ZnCl2/TOTAL_VOLUME*1000.
molarity_NH4Cl = MASS_NH4Cl/MW_NH4Cl/TOTAL_VOLUME*1000.
molarity_KI = MASS_ZnCl2/MW_KI/TOTAL_VOLUME*1000.
molarity_TriEG = MASS_TriEG/MW_TriEG/TOTAL_VOLUME*1000.
molarity_H2O = MASS_H2O/MW_H2O/TOTAL_VOLUME*1000.
if not tqdm_disabled:
print("Electrolyte Composition:")
print("Molarities (moles/L solution): {:.2f} M ZnCl~2~, {:.2f} M NH~4~Cl, {:.2f} M KI, {:.2f} M triethylene glycol, {:.2f} M H~2~O\n".format(molarity_ZnCl2, molarity_NH4Cl, molarity_KI, molarity_TriEG, molarity_H2O))
print("Molalities (moles/kg solution): {:.2f} m ZnCl~2~, {:.2f} m NH~4~Cl, {:.2f} m KI, {:.2f} m triethylene glycol, {:.2f} m H~2~O\n".format(molality_ZnCl2, molality_NH4Cl, molality_KI, molality_TriEG, molality_H2O))
print("Density approx. {:.1f} g/mL".format(density))
filenames = ["22-03-2026-KPS-4.txt"]
all_data = []
for f in filenames:
if len(all_data) == 0:
if "Potentiostat_project" in f:
all_data.append(pd.read_csv(f))
else:
all_data.append(pd.read_csv(f, delimiter="\t"))
else:
df0 = pd.read_csv(f, delimiter="\t")
df0["Elapsed time(s)"] += all_data[-1]["Elapsed time(s)"].iat[-1]
all_data.append(df0)
df = pd.concat(all_data, ignore_index=True)
df["mean_current"] = df["Current(A)"].rolling(2).mean()
df["prev_current"] = df["mean_current"].shift(-1)
df["VChange"] = df["Potential(V)"].diff().abs()
df["is_change"] = (
((df["mean_current"] > 0) & (df["prev_current"] < 0))
| ((df["mean_current"] < 0) & (df["prev_current"] > 0))
).astype(int)
idx_changes = list(df[df["is_change"] == 1].index)
idx_changes.append(len(df) - 1)
all_curves = []
idx_start = 0
for idx in tqdm(idx_changes, disable=tqdm_disabled):
if len(df.iloc[idx_start:idx, :]) > 50:
all_curves.append(df.iloc[idx_start:idx, :])
idx_start = idx
results = []
n_curves = np.max([1, int(np.floor(len(all_curves) / 2))])
for CN in notebook.tnrange(n_curves, disable=tqdm_disabled):
CURVE_N1 = CN * 2
CURVE_N2 = CN * 2 + 1
# Process charge data
if sampling:
N_TERM_POINTS = int(np.min([MIN_POINTS0, len(all_curves[CURVE_N1]) / 2.0]))
MIN_POINTS = int(
np.min([MIN_POINTS0, len(all_curves[CURVE_N1]) - N_TERM_POINTS * 2])
)
df0 = pd.concat(
[
all_curves[CURVE_N1].iloc[:N_TERM_POINTS],
all_curves[CURVE_N1]
.iloc[N_TERM_POINTS:-N_TERM_POINTS]
.sample(n=MIN_POINTS),
all_curves[CURVE_N1].iloc[-N_TERM_POINTS:],
]
).sort_values("Elapsed time(s)", ascending=True)
df0 = df0[df0["VChange"] < DIFF_LIMIT]
else:
df0 = all_curves[CURVE_N1].copy()
df0["mAh"] = np.abs(
scipy.integrate.cumulative_trapezoid(
df0["Current(A)"], df0["Elapsed time(s)"], initial=0
)
* 1000.0
/ 3600.0
)
total_energy0 = scipy.integrate.cumulative_trapezoid(
df0["Current(A)"].abs() * df0["Potential(V)"],
df0["Elapsed time(s)"],
initial=0.0,
)[-1]
# Process discharge data
if sampling:
N_TERM_POINTS = int(np.min([MIN_POINTS0, len(all_curves[CURVE_N2]) / 2.0]))
MIN_POINTS = int(
np.min([MIN_POINTS0, len(all_curves[CURVE_N2]) - N_TERM_POINTS * 2])
)
df1 = pd.concat(
[
all_curves[CURVE_N2].iloc[:N_TERM_POINTS],
all_curves[CURVE_N2]
.iloc[N_TERM_POINTS:-N_TERM_POINTS]
.sample(n=MIN_POINTS),
all_curves[CURVE_N2].iloc[-N_TERM_POINTS:],
]
).sort_values("Elapsed time(s)", ascending=True)
df1 = df1[df1["VChange"] < DIFF_LIMIT]
else:
df1 = all_curves[CURVE_N2].copy()
df1["mAh"] = np.abs(
scipy.integrate.cumulative_trapezoid(
df1["Current(A)"], df1["Elapsed time(s)"], initial=0.0
)
* 1000.0
/ 3600.0
)
total_energy1 = scipy.integrate.cumulative_trapezoid(
df1["Current(A)"].abs() * df1["Potential(V)"],
df1["Elapsed time(s)"],
initial=0.0,
)[-1]
CE = 100.0 * (df1["mAh"].iloc[-1] / df0["mAh"].iloc[-1])
EE = 100.0 * (total_energy1 / total_energy0)
VE = 100.0 * EE / CE
results.append(
{
"Number": CN + 1,
"CE": CE,
"VE": VE,
"EE": EE,
"Charge_potential": df0["Potential(V)"].mean(),
"Discharge_potential": df1["Potential(V)"].mean(),
"Charge_stored": df1["mAh"].iloc[-1] / TOTAL_VOLUME,
"Energy_density_discharge": total_energy1 / TOTAL_VOLUME / 3600.0 * 1000,
}
)
# Save the modified DataFrames back to the all_curves list
all_curves[CURVE_N1] = df0
all_curves[CURVE_N2] = df1
results_df = pd.DataFrame(results)
# results_df = results_df[:4] #only first 4 cycles are worth plotting
if not tqdm_disabled:
print(results_df)
print("")
print(results_df.mean())
# Plot charge/discharge curves
fig1 = go.Figure()
for CN in range(n_curves):
CURVE_N1 = CN * 2
CURVE_N2 = CN * 2 + 1
fig1.add_trace(
go.Scatter(
x=all_curves[CURVE_N1]["mAh"] / TOTAL_VOLUME,
y=all_curves[CURVE_N1]["Potential(V)"],
mode="lines",
name=f"Charge {CN+1}",
line=dict(color="blue", dash="solid"),
)
)
fig1.add_trace(
go.Scatter(
x=all_curves[CURVE_N2]["mAh"] / TOTAL_VOLUME,
y=all_curves[CURVE_N2]["Potential(V)"],
mode="lines",
name=f"Discharge {CN+1}",
line=dict(color="grey", dash="solid"),
)
)
fig1.update_layout(
title="Charge and Discharge Curves",
xaxis_title="Capacity (Ah/L)",
yaxis_title="Potential (V)",
)
fig1.show()
