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
electrolyte = "Cu/Mn"
sampling = True
MIN_POINTS0 = 500
DIFF_LIMIT = 0.1
VOLUME_TO_RESERVOIRS = 3 # unrecorded this time - electrolyte volume in mL, approx, measured by syring when adding to reservoirs
MASS_TO_RESERVOIRS = 3.60 # g of electrolyte actually loaded into system, based on weighing syringe before/after loading reservoirs
density = MASS_TO_RESERVOIRS / VOLUME_TO_RESERVOIRS
# molecular weights in g/mol
# electrolyte component masses, in g
### Cu/Mn ###
MASS_CuSO4 = 2.00 # PENTAHYDRATE
MASS_MnSO4 = 1.29 # MONOHYDRATE
MASS_MSA_99pct = 2.88 # >99.0% w/w solution
MASS_FeSO4 = 0.00 # HEPTAHYDRATE - forgot to add it, should be 0.1 g
MASS_H2O = 7 + 1.57 # 1.5 mL more than Daniel said to get all dissolved
total_mass_kg = (MASS_CuSO4 + MASS_MnSO4 + MASS_MSA_99pct + MASS_FeSO4 + MASS_H2O) / 1000.0
MW_CuSO4 = 249.685 # PENTAHYDRATE
MW_MnSO4 = 169.02 # MONOHYDRATE
MW_MSA = 96.10
MW_FeSO4 = 278.02 # HEPTAHYDRATE
MW_H2O = 18.01528
molality_CuSO4 = MASS_CuSO4 / MW_CuSO4 / total_mass_kg
molality_MnSO4 = MASS_MnSO4 / MW_MnSO4 / total_mass_kg
molality_FeSO4 = MASS_FeSO4 / MW_FeSO4 / total_mass_kg
molality_MSA = MASS_MSA_99pct*0.99 / MW_MSA / total_mass_kg
molality_H2O = (MASS_H2O + 0.01*MASS_MSA_99pct) / MW_H2O / total_mass_kg
molarity_CuSO4 = MASS_CuSO4 / MW_CuSO4 / (total_mass_kg/density)
molarity_MnSO4 = MASS_MnSO4 / MW_MnSO4 / (total_mass_kg/density)
molarity_FeSO4 = MASS_FeSO4 / MW_FeSO4 / (total_mass_kg/density)
molarity_MSA = MASS_MSA_99pct*0.99 / MW_MSA / (total_mass_kg/density)
molarity_H2O = (MASS_H2O + 0.01*MASS_MSA_99pct) / MW_H2O / (total_mass_kg/density)
### ALL-IRON ###
# MASS_FeCl2_2H2O = 1.60
# MASS_MgCl2 = 4.30
# MASS_AscorbicAcid = 0.06 # wasn't able to properly record it
# MASS_HCl_15pct = 0.06 # 15% w/w solution
# MASS_H2O = 9.49
# total_mass_kg = (MASS_FeCl2_2H2O + MASS_MgCl2 + MASS_AscorbicAcid + MASS_HCl_15pct + MASS_H2O) / 1000.0
# MW_FeCl2_2H2O = 162.78
# MW_MgCl2 = 95.211
# MW_AscorbicAcid = 176.124
# MW_HCl = 36.46
# MW_H2O = 18.01528
# molality_FeCl2 = MASS_FeCl2_2H2O / MW_FeCl2_2H2O / total_mass_kg
# molality_MgCl2 = MASS_MgCl2 / MW_MgCl2 / total_mass_kg
# molality_AscorbicAcid = MASS_AscorbicAcid / MW_AscorbicAcid / total_mass_kg
# molality_HCl = MASS_HCl_15pct*0.15 / MW_HCl / total_mass_kg
# molality_H2O = (MASS_H2O + 0.85*MASS_HCl_15pct) / MW_H2O / total_mass_kg
# molarity_FeCl2 = MASS_FeCl2_2H2O / MW_FeCl2_2H2O / VOLUME_TO_RESERVOIRS * 1000.0
# molarity_MgCl2 = MASS_MgCl2 / MW_MgCl2 / VOLUME_TO_RESERVOIRS * 1000.0
# molarity_AscorbicAcid = MASS_AscorbicAcid / MW_AscorbicAcid / VOLUME_TO_RESERVOIRS * 1000.0
# molarity_HCl = MASS_HCl_15pct*0.15 / MW_HCl / VOLUME_TO_RESERVOIRS * 1000.0
# molarity_H2O = (MASS_H2O + 0.85*MASS_HCl_15pct) / MW_H2O / VOLUME_TO_RESERVOIRS * 1000.0
### ZINC-IODIDE ###
# MASS_ZnCl2 = 1.37
# MASS_NH4Cl = 1.06
# MASS_KI = 3.31
# MASS_H2O = 8.52
# MASS_TriEG = 0.63
# total_mass_kg = (MASS_ZnCl2 + MASS_KI + MASS_H2O + MASS_TriEG) / 1000.0
# 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_KI / 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 / VOLUME_TO_RESERVOIRS * 1000.0
# molarity_NH4Cl = MASS_NH4Cl / MW_NH4Cl / VOLUME_TO_RESERVOIRS * 1000.0
# molarity_KI = MASS_KI / MW_KI / VOLUME_TO_RESERVOIRS * 1000.0
# molarity_TriEG = MASS_TriEG / MW_TriEG / VOLUME_TO_RESERVOIRS * 1000.0
# molarity_H2O = MASS_H2O / MW_H2O / VOLUME_TO_RESERVOIRS * 1000.0
filenames = [
"03-06-2026-KPS-28.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:")
if electrolyte=="zinc":
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
)
)
elif electrolyte=="iron":
print(
"Molarities (moles/L solution): {:.2f} M FeCl~2~, {:.2f} M MgCl~2~, {:.2f} M ascorbic acid, {:.2f} M HCl, {:.2f} M H~2~O\n".format(
molarity_FeCl2, molarity_MgCl2, molarity_AscorbicAcid, molarity_HCl, molarity_H2O
)
)
print(
"Molalities (moles/kg solution): {:.2f} m FeCl~2~, {:.2f} m MgCl~2~, {:.2f} m ascorbic acid, {:.2f} m HCl, {:.2f} m H~2~O\n".format(
molality_FeCl2, molarity_MgCl2, molarity_AscorbicAcid, molarity_HCl, molarity_H2O
)
)
elif electrolyte=="Cu/Mn":
print(
"Molarities (moles/L solution): {:.2f} M CuSO~4~, {:.2f} M MnSO~4~, {:.2f} M FeSO~4~, {:.2f} M MSA, {:.2f} M H~2~O\n".format(
molarity_CuSO4, molarity_MnSO4, molarity_FeSO4, molarity_MSA, molarity_H2O
)
)
print(
"Molalities (moles/kg solution): {:.2f} m CuSO~4~, {:.2f} m MnSO~4~, {:.2f} m FeSO~4~, {:.2f} m MSA, {:.2f} m H~2~O\n".format(
molality_CuSO4, molality_MnSO4, molality_FeSO4, molality_MSA, molality_H2O
)
)
else:
print("WARNING: Specify the type of electrolyte")
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(4).mean()
df["prev_current"] = df["mean_current"].shift(-1)
df["VChange"] = df["Potential(V)"].diff().abs()
df["is_OCP"] = (abs(df["mean_current"])<0.0001).astype(int)
df["is_change"] = (
(((df["mean_current"] > 0) & (df["prev_current"] < 0))
| ((df["mean_current"] < 0) & (df["prev_current"] > 0)))
#& df["is_OCP"]==False #current must be greater than 0.5 mA to avoid OCP holds
).astype(int)
df["is_start_charge"] = ((df["prev_current"] < -0.001) & (df["is_OCP"]== 1)).astype(int)
df["is_start_discharge"] = ((df["prev_current"] > 0.001) & (df["is_OCP"]== 1)).astype(int)
idx_changes = list(df[(df["is_start_charge"] == 1) | (df["is_start_discharge"] == 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))])
n_curves = 18
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] / VOLUME_TO_RESERVOIRS,
"Energy_density_discharge": total_energy1 / VOLUME_TO_RESERVOIRS / 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())
