Calculating mass activity and specific activity of Pt/C catalyst on RDE from LSV raw data.
import os
import numpy as np
import pandas as pd
# ============================================================
# 1. CONSTANTS
# ============================================================
DISK_AREA_CM2 = 0.196 # cm^2, change if your RDE disk is different
LOW_POTENTIAL_CUTOFF_V = 0.20
BACKGROUND_POINTS = 10
AVERAGE_WINDOW_HALF_WIDTH = 2 # point itself + 2 below + 2 above = 5-point average
INPUT_FILE = input("Enter the path to the raw LSV text file: ").strip().strip('"').strip("'")
PT_MASS_UG = float(input("Enter Pt mass on disk in micrograms (ug): ").strip())
RESISTANCE_OHM = float(input("Enter the ohmic resistance in ohm: ").strip())
ECSA_M2_PER_G = float(input("Enter the ECSA in m^2/g_Pt: ").strip())
# ============================================================
# 2. OUTPUT FILE NAMES
# ============================================================
input_dir = os.path.dirname(os.path.abspath(INPUT_FILE))
OUTPUT_FILE = os.path.join(input_dir, "processed_lsv.txt")
SUMMARY_FILE = os.path.join(input_dir, "processed_lsv_summary.txt")
# ============================================================
# 3. INPUT CONVERSIONS
# ============================================================
PT_MASS_G = PT_MASS_UG * 1e-6 # ug -> g
if PT_MASS_G <= 0:
raise ValueError("Pt mass must be > 0 g.")
if ECSA_M2_PER_G <= 0:
raise ValueError("ECSA must be > 0 m^2/g.")
if DISK_AREA_CM2 <= 0:
raise ValueError("Disk area must be > 0 cm^2.")
# ============================================================
# 4. READ THE RAW TEXT FILE
# ============================================================
# Assumes the file has a header and two columns:
# first column = potential in V
# second column = current in mA
#
# sep=None with engine='python' lets pandas infer tab/space separators.
df = pd.read_csv(INPUT_FILE, sep=None, engine="python")
if df.shape[1] < 2:
raise ValueError("The input file must contain at least two columns.")
# Keep only the first two columns and rename them clearly
df = df.iloc[:, :2].copy()
df.columns = ["Potential_V_raw", "Current_mA_raw"]
# Convert to numeric in case the parser reads strings
df["Potential_V_raw"] = pd.to_numeric(df["Potential_V_raw"], errors="coerce")
df["Current_mA_raw"] = pd.to_numeric(df["Current_mA_raw"], errors="coerce")
# Drop invalid rows
df = df.dropna(subset=["Potential_V_raw", "Current_mA_raw"]).copy()
if df.empty:
raise ValueError("No valid numeric data found in the file.")
# Sort by potential, because processing steps assume low -> high potential
df = df.sort_values("Potential_V_raw").reset_index(drop=True)
# ============================================================
# 5. CUT NOISY LOW-POTENTIAL DATA
# ============================================================
df = df[df["Potential_V_raw"] >= LOW_POTENTIAL_CUTOFF_V].copy().reset_index(drop=True)
if len(df) < max(BACKGROUND_POINTS, 5):
raise ValueError("Not enough data points remain after cutting below 0.20 V.")
# ============================================================
# 6. APPLY OHMIC RESISTANCE CORRECTION
# ============================================================
# Current is given in mA, so first convert to A:
df["Current_A_raw"] = df["Current_mA_raw"] * 1e-3
# Ohmic loss = I * R, in volts
df["Ohmic_loss_V"] = df["Current_A_raw"] * RESISTANCE_OHM
# iR-corrected potential:
# User requested: corrected potential = raw potential - ohmic loss
df["Potential_V_iR"] = df["Potential_V_raw"] - df["Ohmic_loss_V"]
# ============================================================
# 7. CAPACITIVE BACKGROUND SUBTRACTION
# ============================================================
# Use the last 10 current values after sorting by potential
background_current_mA = df["Current_mA_raw"].tail(BACKGROUND_POINTS).mean()
# Subtract this average from all current values
df["Current_mA_capcorr"] = df["Current_mA_raw"] - background_current_mA
df["Current_A_capcorr"] = df["Current_mA_capcorr"] * 1e-3
# ============================================================
# 8. HELPER FUNCTION FOR LOCAL AVERAGING
# ============================================================
def average_around_target(x_values, y_values, target, half_width=2):
"""
Find the point closest to 'target' in x_values and average y_values
over that point ± half_width.
Returns:
idx_closest, x_closest, y_avg
"""
x_values = np.asarray(x_values)
y_values = np.asarray(y_values)
idx = int(np.argmin(np.abs(x_values - target)))
start = max(0, idx - half_width)
end = min(len(y_values), idx + half_width + 1)
y_avg = np.mean(y_values[start:end])
x_closest = x_values[idx]
return idx, x_closest, y_avg
# ============================================================
# 9. CURRENT AT 0.9 V AND DIFFUSION LIMITING CURRENT AT 0.3 V
# ============================================================
# Use iR-corrected potential to locate the relevant electrochemical potentials
potentials = df["Potential_V_iR"].to_numpy()
currents_A = df["Current_A_capcorr"].to_numpy()
idx_09, pot_09, current_09_A = average_around_target(
potentials, currents_A, target=0.90, half_width=AVERAGE_WINDOW_HALF_WIDTH
)
idx_03, pot_03, current_03_A = average_around_target(
potentials, currents_A, target=0.30, half_width=AVERAGE_WINDOW_HALF_WIDTH
)
# ORR currents are cathodic and usually negative.
# For activity reporting, it is standard to use positive magnitudes.
i_09_abs_A = abs(current_09_A)
i_lim_abs_A = abs(current_03_A)
if i_lim_abs_A <= i_09_abs_A:
raise ValueError(
"The absolute diffusion-limiting current is not larger than the absolute current at 0.9 V. "
"Koutecky-Levich correction would be non-physical with these values."
)
# ============================================================
# 10. KOUTECKY-LEVICH CORRECTION
# ============================================================
# Using magnitudes:
# 1 / i_k = 1 / i - 1 / i_lim
# therefore:
# i_k = (i * i_lim) / (i_lim - i)
kinetic_current_A = (i_09_abs_A * i_lim_abs_A) / (i_lim_abs_A - i_09_abs_A)
# ============================================================
# 11. MASS ACTIVITY
# ============================================================
mass_activity_A_per_g = kinetic_current_A / PT_MASS_G
# ============================================================
# 12. SPECIFIC ACTIVITY
# ============================================================
# ECSA is in m^2/g
# Convert to cm^2/g: 1 m^2 = 1e4 cm^2
ecsa_cm2_per_g = ECSA_M2_PER_G * 1e4
# Specific activity:
# (A/g) / (cm^2/g) = A/cm^2
# Convert to uA/cm^2 by multiplying by 1e6
specific_activity_uA_per_cm2 = (mass_activity_A_per_g / ecsa_cm2_per_g) * 1e6
# ============================================================
# 13. NORMALIZE CURRENT BY DISK AREA FOR OUTPUT
# ============================================================
df["Current_mA_per_cm2"] = df["Current_mA_capcorr"] / DISK_AREA_CM2
# ============================================================
# 14. WRITE STRUCTURED DATA FILE
# ============================================================
with open(OUTPUT_FILE, "w", encoding="utf-8") as f:
f.write("Potential_iR_corrected_V\tCurrent_cap_corrected_mA_per_cm2\n")
for _, row in df.iterrows():
f.write(f"{row['Potential_V_iR']:.8f}\t{row['Current_mA_per_cm2']:.8f}\n")
# ============================================================
# 15. WRITE SUMMARY FILE
# ============================================================
with open(SUMMARY_FILE, "w", encoding="utf-8") as f:
f.write("LSV Processing Summary\n")
f.write("======================\n\n")
f.write(f"Number_of_data_points_after_cutoff\t{len(df)}\n")
f.write(f"ECSA_m2_per_g\t{ECSA_M2_PER_G:.6f}\n")
f.write(f"Mass_activity_A_per_g\t{mass_activity_A_per_g:.6f}\n")
f.write(f"Specific_activity_uA_per_cm2\t{specific_activity_uA_per_cm2:.6f}\n")
f.write(f"Pt_mass_ug\t{PT_MASS_UG:.6f}\n")
f.write(f"Resistance_ohm\t{RESISTANCE_OHM:.6f}\n")
f.write(f"Disk_area_cm2\t{DISK_AREA_CM2:.6f}\n")
f.write(f"Background_current_mA\t{background_current_mA:.6f}\n")
f.write(f"Closest_potential_for_0p9V_V\t{pot_09:.6f}\n")
f.write(f"Averaged_current_at_0p9V_A\t{current_09_A:.9e}\n")
f.write(f"Closest_potential_for_0p3V_V\t{pot_03:.6f}\n")
f.write(f"Averaged_diffusion_limiting_current_A\t{current_03_A:.9e}\n")
f.write(f"Kinetic_current_A\t{kinetic_current_A:.9e}\n")
# ============================================================
# 16. PRINT SUMMARY TO SCREEN
# ============================================================
print("\nProcessing completed successfully.\n")
print(f"Number of data points after cutoff: {len(df)}")
print(f"Background current (last {BACKGROUND_POINTS} points): {background_current_mA:.6f} mA")
print(f"Closest iR-corrected potential to 0.90 V: {pot_09:.6f} V")
print(f"Averaged current at ~0.90 V: {current_09_A:.9e} A")
print(f"Closest iR-corrected potential to 0.30 V: {pot_03:.6f} V")
print(f"Averaged diffusion-limiting current at ~0.30 V: {current_03_A:.9e} A")
print(f"Kinetic current at 0.90 V: {kinetic_current_A:.9e} A")
print(f"Mass activity: {mass_activity_A_per_g:.6f} A/g")
print(f"Specific activity: {specific_activity_uA_per_cm2:.6f} uA/cm^2")
print(f"\nStructured data written to: {OUTPUT_FILE}")
print(f"Summary written to: {SUMMARY_FILE}")