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
from plotly.subplots import make_subplots
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.38
MASS_NH4Cl = 1.06
MASS_KI = 3.32
MASS_H2O = 8.50
MASS_TriEG = 0.55
total_mass_kg = (MASS_ZnCl2 + MASS_KI + MASS_H2O + MASS_TriEG) / 1000.0
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.60 # 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.0
molarity_NH4Cl = MASS_NH4Cl / MW_NH4Cl / TOTAL_VOLUME * 1000.0
molarity_KI = MASS_ZnCl2 / MW_KI / TOTAL_VOLUME * 1000.0
molarity_TriEG = MASS_TriEG / MW_TriEG / TOTAL_VOLUME * 1000.0
molarity_H2O = MASS_H2O / MW_H2O / TOTAL_VOLUME * 1000.0
filenames = [
"16-04-2026-KPS-22.zip",
]
all_data = []
for f in filenames:
if len(all_data) == 0:
all_data.append(pd.read_csv(f, delimiter="\t").dropna())
else:
df0 = pd.read_csv(f, delimiter="\t").dropna()
df0["Elapsed time(s)"] += all_data[-1]["Elapsed time(s)"].iat[-1]
all_data.append(df0)
df = pd.concat(all_data, ignore_index=True)
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\n".format(density))
print(
"Experiment length: {:.1f} hours".format(
all_data[-1]["Elapsed time(s)"].iat[-1] / 3600.0
)
)
df["mean_current"] = df["Current(A)"].rolling(3).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)
if not tqdm_disabled:
print(results_df)
print("")
print(results_df.mean())
# Color gradient for charge curves
charge_colors = [
px.colors.sequential.Blues[int(i)]
for i in np.linspace(3, len(px.colors.sequential.Blues) - 1, n_curves)
]
discharge_colors = [
px.colors.sequential.Greys[int(i)]
for i in np.linspace(3, len(px.colors.sequential.Greys) - 1, n_curves)
]
# 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=charge_colors[CN], dash="solid"),
showlegend=False,
)
)
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=discharge_colors[CN], dash="solid"),
showlegend=False,
)
)
fig1.update_layout(
xaxis_title="Capacity (Ah/L)",
yaxis_title="Potential (V)",
legend=dict(orientation="h", yanchor="bottom", y=0.02, xanchor="right", x=0.99),
hoverlabel=dict(
bgcolor="white",
),
xaxis=dict(range=[-1, 10]),
yaxis=dict(range=[-.49, 1.8]),
)
fig1.add_trace(
go.Scatter(
x=[None],
y=[None],
mode="lines",
line=dict(color=charge_colors[0], dash="solid"),
name="Charge (cycle 1)",
)
)
fig1.add_trace(
go.Scatter(
x=[None],
y=[None],
mode="lines",
line=dict(color=charge_colors[-1], dash="solid"),
name=f"Charge (cycle {n_curves})",
)
)
fig1.add_trace(
go.Scatter(
x=[None],
y=[None],
mode="lines",
line=dict(color=discharge_colors[0], dash="solid"),
name="Discharge (cycle 1)",
)
)
fig1.add_trace(
go.Scatter(
x=[None],
y=[None],
mode="lines",
line=dict(color=discharge_colors[-1], dash="solid"),
name=f"Discharge (cycle {n_curves})",
)
)
fig1.show()
