ECSA_Linear_Background_Hupd

import os

# ============================================================
# 1. CONSTANTS
# ============================================================
SCAN_RATE_MV_S = 20.0
SCAN_RATE_V_S = SCAN_RATE_MV_S / 1000.0   # V/s
H_CHARGE = 2.1e-4                         # C/cm^2_Pt
TARGET_MAX_POTENTIAL = 1.23               # V
MIN_POTENTIAL_WINDOW_LOW = 0.04
MIN_POTENTIAL_WINDOW_HIGH = 0.06


# ============================================================
# 2. USER INPUTS
# ============================================================
R_ohm = float(input("Enter the ohmic resistance R (ohm): "))
mass_ug = float(input("Enter the Pt mass on the electrode (ug): "))
mass_g = mass_ug * 1e-6

file_path = input("Please enter the full path to the raw CV text file: ").strip()
file_path = file_path.strip('"').strip("'")


# ============================================================
# 3. READ RAW CV FILE
# ============================================================
potential = []
current_mA = []

directory = os.path.dirname(file_path)
filename = os.path.basename(file_path)
name, extension = os.path.splitext(filename)

try:
    with open(file_path, "r") as f:
        lines = f.readlines()

    for line in lines:
        if "Ewe/V" in line:
            continue

        line = line.replace("\t", " ").strip()
        if not line:
            continue

        parts = line.split()
        if len(parts) < 2:
            continue

        potential.append(float(parts[0]))
        current_mA.append(float(parts[1]))

except FileNotFoundError:
    raise FileNotFoundError("The specified file was not found. Please check the file path and try again.")
except Exception as e:
    raise RuntimeError(f"An error occurred while reading the input file: {e}")

if len(potential) < 10:
    raise ValueError("Not enough data points were found in the file.")


# ============================================================
# 4. OPTIONAL POTENTIAL RANGE CORRECTION
# ============================================================
potential_correction = max(potential) - TARGET_MAX_POTENTIAL
potential = [p - potential_correction for p in potential]


# ============================================================
# 5. FIND THE START OF THE THIRD CV CYCLE
# ============================================================
cycle_count = 0
start_index_last_cycle = None

for i in range(2, len(potential) - 2):
    if i > 150:
        is_local_min = (
            potential[i] < potential[i - 1]
            and potential[i] < potential[i + 1]
            and potential[i] < potential[i + 2]
        )
        is_in_low_window = (MIN_POTENTIAL_WINDOW_LOW < potential[i] < MIN_POTENTIAL_WINDOW_HIGH)

        if is_local_min and is_in_low_window:
            cycle_count += 1

        if cycle_count == 2:
            start_index_last_cycle = i
            break

if start_index_last_cycle is None:
    raise ValueError("Could not identify the start of the third CV cycle.")


# ============================================================
# 6. EXTRACT LAST CYCLE
# ============================================================
potential_last = potential[start_index_last_cycle:]
current_last_mA = current_mA[start_index_last_cycle:]

if len(potential_last) < 10:
    raise ValueError("The extracted last cycle is too short.")


# ============================================================
# 7. SMOOTH EDGES / REMOVE VERY ABNORMAL POINTS
# ============================================================
potential_last_smoothed = []
current_last_smoothed_mA = []

for i in range(len(potential_last)):
    if i == 0:
        potential_last_smoothed.append(potential_last[i])
        current_last_smoothed_mA.append(current_last_mA[i])

    elif 0 < i < len(potential_last) - 1:
        prev_abs = abs(current_last_mA[i - 1])
        next_abs = abs(current_last_mA[i + 1])

        if prev_abs == 0 or next_abs == 0:
            potential_last_smoothed.append(potential_last[i])
            current_last_smoothed_mA.append(current_last_mA[i])
        else:
            var_prev = abs(current_last_mA[i] - current_last_mA[i - 1]) / prev_abs
            var_next = abs(current_last_mA[i] - current_last_mA[i + 1]) / next_abs

            if var_prev < 0.3 and var_next < 0.3:
                potential_last_smoothed.append(potential_last[i])
                current_last_smoothed_mA.append(current_last_mA[i])

# close the cycle visually
potential_last_smoothed.append(potential_last_smoothed[0])
current_last_smoothed_mA.append(current_last_smoothed_mA[0])

if len(potential_last_smoothed) < 10:
    raise ValueError("Too few points remain after smoothing/filtering.")


# ============================================================
# 8. APPLY iR CORRECTION TO LAST CYCLE
# ============================================================
potential_ir_corrected = []
for p, i_mA in zip(potential_last_smoothed, current_last_smoothed_mA):
    p_corr = p - R_ohm * i_mA / 1000.0
    potential_ir_corrected.append(p_corr)


# ============================================================
# 9. COMPUTE ECSA FROM iR-CORRECTED LAST CYCLE
# ============================================================
index_right = None
for i in range(len(potential_ir_corrected)):
    if potential_ir_corrected[i] >= 0.45:
        index_right = i
        break

if index_right is None:
    raise ValueError("Could not find the right integration limit at 0.45 V in the iR-corrected cycle.")

index_left = None
current_right = current_last_smoothed_mA[index_right]

for i in range(len(current_last_smoothed_mA)):
    if current_last_smoothed_mA[i] >= current_right:
        index_left = i
        break

if index_left is None or index_left >= index_right:
    raise ValueError("Could not identify a valid Hupd integration window.")

potential_to_integrate = potential_ir_corrected[index_left:index_right]
current_to_integrate_mA = current_last_smoothed_mA[index_left:index_right]

if len(potential_to_integrate) < 2:
    raise ValueError("Not enough points in the integration range.")

integrated_area = 0.0
for i in range(len(potential_to_integrate) - 1):
    dE = potential_to_integrate[i + 1] - potential_to_integrate[i]
    integrated_area += (current_to_integrate_mA[i] / 1000.0) * dE

capacitive_background = (
    current_to_integrate_mA[-1] / 1000.0
) * (potential_to_integrate[-1] - potential_to_integrate[0])

integrated_area -= capacitive_background

ecsa_m2_g = integrated_area / (SCAN_RATE_V_S * H_CHARGE * mass_g * 10000.0)


# ============================================================
# 10. SAVE OUTPUT FILES
# ============================================================
data_filename = f"{name}_iR-corrected_last-cycle{extension}"
data_file_path = os.path.join(directory, data_filename)

summary_filename = f"{name}_ECSA_summary.txt"
summary_file_path = os.path.join(directory, summary_filename)

try:
    # Clean numeric data file
    with open(data_file_path, "w") as f_data:
        f_data.write("Ewe/V\t<I>/mA\n")
        for p, i_mA in zip(potential_ir_corrected, current_last_smoothed_mA):
            f_data.write(f"{p:.8f}\t{i_mA:.8f}\n")

    # Separate summary file
    with open(summary_file_path, "w") as f_summary:
        f_summary.write("ECSA ANALYSIS SUMMARY\n")
        f_summary.write("=====================\n")
        f_summary.write(f"Input raw CV file: {file_path}\n")
        f_summary.write(f"Corrected CV output file: {data_file_path}\n")
        f_summary.write(f"Ohmic resistance (ohm): {R_ohm:.6f}\n")
        f_summary.write(f"Pt mass on electrode (ug): {mass_ug:.6f}\n")
        f_summary.write(f"Pt mass on electrode (g): {mass_g:.10e}\n")
        f_summary.write(f"Scan rate (mV/s): {SCAN_RATE_MV_S:.3f}\n")
        f_summary.write(f"Scan rate (V/s): {SCAN_RATE_V_S:.6f}\n")
        f_summary.write(f"Hupd charge density (C/cm^2_Pt): {H_CHARGE:.6e}\n")
        f_summary.write(f"Integration start index: {index_left}\n")
        f_summary.write(f"Integration end index: {index_right}\n")
        f_summary.write(f"Integration start potential (V): {potential_ir_corrected[index_left]:.8f}\n")
        f_summary.write(f"Integration end potential (V): {potential_ir_corrected[index_right]:.8f}\n")
        f_summary.write(f"Integrated charge after background subtraction: {integrated_area:.10e}\n")
        f_summary.write(f"ECSA (m^2/g_Pt): {ecsa_m2_g:.6f}\n")

    print(f"Corrected CV data file created: {data_file_path}")
    print(f"Summary file created: {summary_file_path}")
    print(f"ECSA = {ecsa_m2_g:.4f} m^2/g_Pt")

except Exception as e:
    raise RuntimeError(f"An error occurred while creating the output files: {e}")