ECSA_CO+Hupd_Final

# ============================================================
# 1) CONSTANTS
# ============================================================
import numpy as np
import pandas as pd
from pathlib import Path

# Expected specific charges
# CO monolayer oxidation on Pt: 420 microC cm^-2_Pt
# Hupd on polycrystalline Pt:   210 microC cm^-2_Pt
CO_STRIPPING_CHARGE_DENSITY_C_PER_CM2 = 420e-6
HUPD_CHARGE_DENSITY_C_PER_CM2 = 210e-6

# File parsing
DEFAULT_SKIP_INITIAL_ROWS = 3
POTENTIAL_COLUMN_CANDIDATES = ["Ewe/V", "Potential_V", "Potential", "E/V"]
CURRENT_COLUMN_CANDIDATES = ["<I>/mA", "Current_mA", "Current", "I/mA"]

# Cycle detection
LOW_POTENTIAL_TARGET_V = 0.05
LOW_POTENTIAL_TOLERANCE_V = 0.08
MIN_NEGATIVE_STEPS = 2
MIN_POSITIVE_STEPS = 2
MIN_CYCLE_SEPARATION_ROWS = 200

# CO stripping integration
INTERSECTION_SEARCH_MIN_V = 0.50
CO_RIGHT_LIMIT_TARGET_V = 1.22

# Hupd integration
HUPD_RIGHT_LIMIT_TARGET_V = 0.45

# Electrochemistry constants
SCAN_RATE_MV_S = 20.0
SCAN_RATE_V_S = SCAN_RATE_MV_S / 1000.0

# Output formatting
OUTPUT_FLOAT_FORMAT = "%.8f"

# ============================================================
# 2) INPUT VALUES
# ============================================================
file_path = input("Enter the full path to the raw CV text file: ").strip().strip('"').strip("'")
r_ohm = float(input("Enter the ohmic resistance in ohms for iR correction: ").strip())

pt_mass_ug = float(input("Enter the Pt mass on the electrode in micrograms (ug): ").strip())
pt_mass_g = pt_mass_ug * 1e-6

if pt_mass_g <= 0:
    raise ValueError("The Pt mass must be positive.")

skip_initial_rows = DEFAULT_SKIP_INITIAL_ROWS

# ============================================================
# 3) TEXT FILES TO ANALYZE
# ============================================================
input_file = Path(file_path)
if not input_file.exists():
    raise FileNotFoundError(f"File not found: {input_file}")

def find_existing_column(columns, candidates):
    for cand in candidates:
        if cand in columns:
            return cand
    raise KeyError(f"Could not find any of these columns: {candidates}. Found columns: {list(columns)}")

df_raw = pd.read_csv(input_file, sep=r"\s+|\t+", engine="python")

potential_col = find_existing_column(df_raw.columns, POTENTIAL_COLUMN_CANDIDATES)
current_col = find_existing_column(df_raw.columns, CURRENT_COLUMN_CANDIDATES)

df = df_raw[[potential_col, current_col]].copy()
df.columns = ["Potential_V_measured", "Current_mA"]
df = df.apply(pd.to_numeric, errors="coerce").dropna().reset_index(drop=True)

if len(df) <= skip_initial_rows + 20:
    raise ValueError("Not enough valid rows after applying the initial row skip.")

df = df.iloc[skip_initial_rows:].reset_index(drop=True)

# ============================================================
# 4) CALCULATIONS
# ============================================================
def moving_average(arr, window=3):
    if window <= 1:
        return np.asarray(arr, dtype=float).copy()
    return pd.Series(arr).rolling(window=window, center=True, min_periods=1).mean().to_numpy()

def detect_cycle_start_indices(potential_array):
    """
    Detect cycle starts from the point where the potential reaches the low-potential
    turning point and starts increasing again.
    """
    p = moving_average(np.asarray(potential_array, dtype=float), window=3)
    dp = np.diff(p)

    starts = [0]
    i = MIN_NEGATIVE_STEPS

    while i < len(dp) - MIN_POSITIVE_STEPS:
        neg_ok = np.all(dp[i - MIN_NEGATIVE_STEPS:i] < 0)
        pos_ok = np.all(dp[i:i + MIN_POSITIVE_STEPS] > 0)

        if neg_ok and pos_ok and p[i] <= LOW_POTENTIAL_TARGET_V + LOW_POTENTIAL_TOLERANCE_V:
            new_start = i
            if new_start - starts[-1] > MIN_CYCLE_SEPARATION_ROWS:
                starts.append(new_start)
            i += 20
        else:
            i += 1

    return starts

cycle_starts = detect_cycle_start_indices(df["Potential_V_measured"].to_numpy())

if len(cycle_starts) < 3:
    raise ValueError(
        f"Only {len(cycle_starts)} cycle starts were detected. "
        "Please check the file format or adjust the cycle detection settings."
    )

cycle_starts = cycle_starts[:3]
cycle_bounds = []
for idx, start in enumerate(cycle_starts):
    end = cycle_starts[idx + 1] if idx + 1 < len(cycle_starts) else len(df)
    cycle_bounds.append((start, end))

cycles = []
for start, end in cycle_bounds:
    cycle_df = df.iloc[start:end].reset_index(drop=True)
    cycles.append(cycle_df)

if len(cycles) < 3:
    raise ValueError("Three cycles were expected but could not be isolated.")

cycle_1 = cycles[0].copy()  # contains CO stripping peak
cycle_2 = cycles[1].copy()
cycle_3 = cycles[2].copy()  # CO background; also used for Hupd signal

def split_anodic_branch(cycle_df):
    """Keep only the forward scan from low potential up to the maximum potential."""
    idx_max = int(cycle_df["Potential_V_measured"].idxmax())
    return cycle_df.iloc[:idx_max + 1].reset_index(drop=True)

cycle_1_anodic = split_anodic_branch(cycle_1)
cycle_3_anodic = split_anodic_branch(cycle_3)

# iR correction: E_corrected = E_measured - R * I
# Current in the file is in mA, so convert to A before multiplying by ohms
for c in (cycle_1, cycle_2, cycle_3, cycle_1_anodic, cycle_3_anodic):
    current_a = c["Current_mA"].to_numpy() / 1000.0
    c["Potential_V_iR"] = c["Potential_V_measured"].to_numpy() - r_ohm * current_a

def choose_nearest_index_and_value(x, target):
    idx = int(np.argmin(np.abs(np.asarray(x, dtype=float) - float(target))))
    return idx, float(np.asarray(x, dtype=float)[idx])

def trapezoid_integral(x, y):
    """
    Compatibility trapezoidal integral that does not rely on np.trapz.
    Returns the area under y(x).
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)

    if len(x) != len(y):
        raise ValueError("x and y must have the same length for integration.")
    if len(x) < 2:
        return 0.0

    dx = np.diff(x)
    avg_y = 0.5 * (y[:-1] + y[1:])
    return float(np.sum(dx * avg_y))

def integrate_between_limits_on_xy(x_data, y_data, left_v, right_v):
    """
    Integrate y vs x between left_v and right_v using linear interpolation at boundaries.
    """
    x = np.asarray(x_data, dtype=float)
    y = np.asarray(y_data, dtype=float)

    if len(x) != len(y):
        raise ValueError("x and y must have the same length.")
    if len(x) < 2:
        return 0.0, np.array([]), np.array([])

    if left_v < x.min() or right_v > x.max():
        raise ValueError("Integration limits are outside the available potential window.")
    if right_v <= left_v:
        raise ValueError("The right integration limit must be larger than the left limit.")

    x_mid = x[(x > left_v) & (x < right_v)]
    x_int = np.concatenate(([left_v], x_mid, [right_v]))
    y_int = np.interp(x_int, x, y)

    area = trapezoid_integral(x_int, y_int)
    return float(area), x_int, y_int

def find_first_intersection_above_threshold(cycle1_anodic_df, cycle3_anodic_df, min_potential_v=0.50):
    """
    Find the first intersection above the threshold potential by interpolating the
    background current of cycle 3 onto the potential grid of cycle 1.
    Uses iR-corrected potential for the intersection.
    """
    x1 = cycle1_anodic_df["Potential_V_iR"].to_numpy()
    y1 = cycle1_anodic_df["Current_mA"].to_numpy()

    x3 = cycle3_anodic_df["Potential_V_iR"].to_numpy()
    y3 = cycle3_anodic_df["Current_mA"].to_numpy()

    common_min = max(min_potential_v, float(np.min(x1)), float(np.min(x3)))
    common_max = min(float(np.max(x1)), float(np.max(x3)))

    mask = (x1 >= common_min) & (x1 <= common_max)
    x = x1[mask]
    y1_use = y1[mask]

    if len(x) < 2:
        raise ValueError("Not enough overlapping points were found to determine the CO intersection.")

    y3_interp = np.interp(x, x3, y3)
    delta = y1_use - y3_interp

    for i in range(len(delta) - 1):
        if delta[i] == 0:
            return float(x[i])
        if delta[i] * delta[i + 1] < 0:
            x_left, x_right = x[i], x[i + 1]
            d_left, d_right = delta[i], delta[i + 1]
            x_cross = x_left - d_left * (x_right - x_left) / (d_right - d_left)
            return float(x_cross)

    raise ValueError("No intersection between cycle 1 and cycle 3 was found above the chosen threshold.")

# ---------------- CO stripping ECSA ----------------
co_left_limit_v = find_first_intersection_above_threshold(
    cycle_1_anodic, cycle_3_anodic, min_potential_v=INTERSECTION_SEARCH_MIN_V
)

_, co_right_limit_v = choose_nearest_index_and_value(
    cycle_1_anodic["Potential_V_iR"].to_numpy(), CO_RIGHT_LIMIT_TARGET_V
)

if co_right_limit_v <= co_left_limit_v:
    raise ValueError(
        f"CO right integration limit ({co_right_limit_v:.4f} V) is not larger than the "
        f"left integration limit ({co_left_limit_v:.4f} V)."
    )

area_cycle1_co_mA_V, _, _ = integrate_between_limits_on_xy(
    cycle_1_anodic["Potential_V_iR"].to_numpy(),
    cycle_1_anodic["Current_mA"].to_numpy(),
    co_left_limit_v,
    co_right_limit_v,
)

area_cycle3_co_mA_V, _, _ = integrate_between_limits_on_xy(
    cycle_3_anodic["Potential_V_iR"].to_numpy(),
    cycle_3_anodic["Current_mA"].to_numpy(),
    co_left_limit_v,
    co_right_limit_v,
)

co_peak_area_mA_V = area_cycle1_co_mA_V - area_cycle3_co_mA_V

if SCAN_RATE_V_S <= 0:
    raise ValueError("The scan rate must be positive.")

# mA·V / (V/s) = mC ; then convert mC to C
co_stripping_charge_C = (co_peak_area_mA_V / SCAN_RATE_V_S) / 1000.0

# Geometric Pt area from CO stripping
ecsa_co_cm2_pt = co_stripping_charge_C / CO_STRIPPING_CHARGE_DENSITY_C_PER_CM2

# Mass-normalized ECSA: cm^2/g -> m^2/g by multiplying by 1e-4
ecsa_co_m2_g_pt = (ecsa_co_cm2_pt / pt_mass_g) * 1e-4

# ---------------- Hupd ECSA ----------------
def first_positive_current_potential(cycle_df):
    x = cycle_df["Potential_V_iR"].to_numpy()
    y = cycle_df["Current_mA"].to_numpy()
    pos_idx = np.where(y > 0)[0]
    if len(pos_idx) == 0:
        raise ValueError("No positive-current point was found in the selected anodic branch.")
    idx = int(pos_idx[0])
    return float(x[idx])

hupd_left_cycle1_v = first_positive_current_potential(cycle_1_anodic)
hupd_left_cycle3_v = first_positive_current_potential(cycle_3_anodic)
hupd_left_limit_v = max(hupd_left_cycle1_v, hupd_left_cycle3_v)

_, hupd_right_limit_v = choose_nearest_index_and_value(
    cycle_3_anodic["Potential_V_iR"].to_numpy(), HUPD_RIGHT_LIMIT_TARGET_V
)

if hupd_right_limit_v <= hupd_left_limit_v:
    raise ValueError(
        f"Hupd right integration limit ({hupd_right_limit_v:.4f} V) is not larger than the "
        f"left integration limit ({hupd_left_limit_v:.4f} V)."
    )

area_cycle1_hupd_mA_V, _, _ = integrate_between_limits_on_xy(
    cycle_1_anodic["Potential_V_iR"].to_numpy(),
    cycle_1_anodic["Current_mA"].to_numpy(),
    hupd_left_limit_v,
    hupd_right_limit_v,
)

area_cycle3_hupd_mA_V, _, _ = integrate_between_limits_on_xy(
    cycle_3_anodic["Potential_V_iR"].to_numpy(),
    cycle_3_anodic["Current_mA"].to_numpy(),
    hupd_left_limit_v,
    hupd_right_limit_v,
)

# Hupd signal = cycle 3 - cycle 1 background
hupd_net_area_mA_V = area_cycle3_hupd_mA_V - area_cycle1_hupd_mA_V
hupd_charge_C = (hupd_net_area_mA_V / SCAN_RATE_V_S) / 1000.0

ecsa_hupd_cm2_pt = hupd_charge_C / HUPD_CHARGE_DENSITY_C_PER_CM2
ecsa_hupd_m2_g_pt = (ecsa_hupd_cm2_pt / pt_mass_g) * 1e-4

# ============================================================
# 5) OUTPUTS
# ============================================================
output_dir = input_file.parent

cycle1_output = cycle_1[["Potential_V_iR", "Current_mA"]].copy()
cycle3_output = cycle_3[["Potential_V_iR", "Current_mA"]].copy()

cycle1_output.columns = ["Potential_V_iR", "Current_mA"]
cycle3_output.columns = ["Potential_V_iR", "Current_mA"]

cycle1_output_path = output_dir / f"{input_file.stem}_cycle1_CO_stripping_iR_corrected.txt"
cycle3_output_path = output_dir / f"