Pt-C Activity RDE

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 TARGET POTENTIALS 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()

TARGET_ACTIVITY_POTENTIALS_V = [0.90, 0.95]
DIFFUSION_LIMITING_POTENTIAL_V = 0.30

idx_03, pot_03, current_03_A = average_around_target(
    potentials,
    currents_A,
    target=DIFFUSION_LIMITING_POTENTIAL_V,
    half_width=AVERAGE_WINDOW_HALF_WIDTH
)

# ORR currents are cathodic and usually negative.
# For activity reporting, it is standard to use positive magnitudes.
i_lim_abs_A = abs(current_03_A)

activity_results = {}

for target_potential in TARGET_ACTIVITY_POTENTIALS_V:
    idx_target, pot_target, current_target_A = average_around_target(
        potentials,
        currents_A,
        target=target_potential,
        half_width=AVERAGE_WINDOW_HALF_WIDTH
    )

    i_target_abs_A = abs(current_target_A)

    if i_lim_abs_A <= i_target_abs_A:
        raise ValueError(
            f"The absolute diffusion-limiting current is not larger than the absolute current at "
            f"{target_potential:.2f} 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_target_abs_A * i_lim_abs_A) / (i_lim_abs_A - i_target_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

    label = f"{target_potential:.2f}V".replace(".", "p")

    activity_results[label] = {
        "target_potential_V": target_potential,
        "closest_potential_V": pot_target,
        "averaged_current_A": current_target_A,
        "kinetic_current_A": kinetic_current_A,
        "mass_activity_A_per_g": mass_activity_A_per_g,
        "specific_activity_uA_per_cm2": specific_activity_uA_per_cm2,
    }


# ============================================================
# 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"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_0p3V_V\t{pot_03:.6f}\n")
    f.write(f"Averaged_diffusion_limiting_current_A\t{current_03_A:.9e}\n\n")

    f.write("Activity Results\n")
    f.write("================\n\n")

    for label, result in activity_results.items():
        f.write(f"Target_potential_{label}_V\t{result['target_potential_V']:.6f}\n")
        f.write(f"Closest_potential_for_{label}_V\t{result['closest_potential_V']:.6f}\n")
        f.write(f"Averaged_current_at_{label}_A\t{result['averaged_current_A']:.9e}\n")
        f.write(f"Kinetic_current_at_{label}_A\t{result['kinetic_current_A']:.9e}\n")
        f.write(f"Mass_activity_at_{label}_A_per_g\t{result['mass_activity_A_per_g']:.6f}\n")
        f.write(f"Specific_activity_at_{label}_uA_per_cm2\t{result['specific_activity_uA_per_cm2']:.6f}\n\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.30 V: {pot_03:.6f} V")
print(f"Averaged diffusion-limiting current at ~0.30 V: {current_03_A:.9e} A")

for label, result in activity_results.items():
    print("")
    print(f"Activity at {result['target_potential_V']:.2f} V vs RHE")
    print("-" * 40)
    print(f"Closest iR-corrected potential: {result['closest_potential_V']:.6f} V")
    print(f"Averaged current: {result['averaged_current_A']:.9e} A")
    print(f"Kinetic current: {result['kinetic_current_A']:.9e} A")
    print(f"Mass activity: {result['mass_activity_A_per_g']:.6f} A/g")
    print(f"Specific activity: {result['specific_activity_uA_per_cm2']:.6f} uA/cm^2")

print(f"\nStructured data written to: {OUTPUT_FILE}")
print(f"Summary written to: {SUMMARY_FILE}")