Calculates ECSA of Pt from Hupd region assuming linear background, from a text file containing three CV cycles. The last one is identified, iR corrected and used for calculations
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}")