Pt size distribution TEM - NEW

# GetHug creator: Alessio
# Known constants
units = {}

MIN_DIAMETER_NM = 0.00
units["MIN_DIAMETER_NM"] = "nm"

MAX_DIAMETER_NM = 50.00
units["MAX_DIAMETER_NM"] = "nm"

DEFAULT_BIN_WIDTH_NM = 0.25
units["DEFAULT_BIN_WIDTH_NM"] = "nm"

PT_DENSITY_G_PER_CM3 = 21.45
units["PT_DENSITY_G_PER_CM3"] = "g/cm^3"

NM_TO_CM = 1e-7
units["NM_TO_CM"] = "cm/nm"

NM2_TO_M2 = 1e-18
units["NM2_TO_M2"] = "m^2/nm^2"

OUTPUT_FILE = "particle_size_distribution.txt"
units["OUTPUT_FILE"] = "text"

# Correction factor for cuboctahedron
SURFACE_AREA_CORRECTION_FACTOR = 1.105

# Input values
try:
    user_bin_width_text = input(f"Enter bin width in nm (press Enter to use {DEFAULT_BIN_WIDTH_NM} nm): ").strip()
    if user_bin_width_text == "":
        bin_width_nm = DEFAULT_BIN_WIDTH_NM
    else:
        bin_width_nm = float(user_bin_width_text)
    units["bin_width_nm"] = "nm"

    if bin_width_nm <= 0:
        raise ValueError("Bin width must be greater than 0.")

except ValueError as exc:
    raise ValueError(f"Invalid bin width input. Please enter a positive number in nm. Details: {exc}")

# Files to analyze
INPUT_FILE_1 = "file_1.txt"

# Calculations
import math
import statistics

diameters_nm = []

with open(INPUT_FILE_1, "r", encoding="utf-8") as f:
    lines = f.readlines()

header_found = False
length_col_index = None

for line in lines:
    stripped_line = line.strip()

    if not stripped_line:
        continue

    if stripped_line.startswith("#"):
        continue

    parts = stripped_line.split("\t")

    if not header_found:
        normalized_parts = [p.strip() for p in parts]
        if len(normalized_parts) == 6:
            header_found = True
            length_col_index = 6  # Assuming the diameter is in the 7th column (index 6)
        continue

    if length_col_index is None:
        continue

    if len(parts) <= length_col_index:
        continue

    cleaned_parts = [p.strip() for p in parts]

    if any(p == "" for p in cleaned_parts):
        continue

    try:
        diameter_nm = float(cleaned_parts[length_col_index])
        units["diameter_nm"] = "nm"
    except ValueError:
        continue

    if math.isnan(diameter_nm):
        continue

    if diameter_nm < MIN_DIAMETER_NM or diameter_nm > MAX_DIAMETER_NM:
        continue

    diameters_nm.append(diameter_nm)

particle_count = len(diameters_nm)
units["particle_count"] = "count"

if particle_count == 0:
    raise ValueError("No valid particle diameters were found in the file within the selected diameter range.")

number_of_bins = int(math.ceil((MAX_DIAMETER_NM - MIN_DIAMETER_NM) / bin_width_nm))
units["number_of_bins"] = "count"

total_range_nm = MAX_DIAMETER_NM - MIN_DIAMETER_NM
units["total_range_nm"] = "nm"

bin_counts_total = 0
units["bin_counts_total"] = "count"

total_surface_area_m2 = 0.0
units["total_surface_area_m2"] = "m^2"

total_mass_g = 0.0
units["total_mass_g"] = "g"

distribution_rows = []

for bin_index in range(number_of_bins):
    bin_lower_nm = MIN_DIAMETER_NM + bin_index * bin_width_nm
    units["bin_lower_nm"] = "nm"

    bin_upper_nm = bin_lower_nm + bin_width_nm
    units["bin_upper_nm"] = "nm"

    bin_center_nm = bin_lower_nm + (bin_width_nm / 2.0)
    units["bin_center_nm"] = "nm"

    if bin_index == number_of_bins - 1:
        count_in_bin = sum(1 for d in diameters_nm if bin_lower_nm <= d <= MAX_DIAMETER_NM)
    else:
        count_in_bin = sum(1 for d in diameters_nm if bin_lower_nm <= d < bin_upper_nm)
    units["count_in_bin"] = "count"

    frequency_in_bin = count_in_bin / particle_count
    units["frequency_in_bin"] = "fraction"

    bin_counts_total += count_in_bin

    distribution_rows.append((bin_center_nm, count_in_bin, frequency_in_bin))

for diameter_nm_single in diameters_nm:
    radius_nm = diameter_nm_single / 2.0
    units["radius_nm"] = "nm"

    surface_area_nm2 = 4.0 * math.pi * (radius_nm ** 2)
    units["surface_area_nm2"] = "nm^2"

    radius_cm = radius_nm * NM_TO_CM
    units["radius_cm"] = "cm"

    volume_cm3 = (4.0 / 3.0) * math.pi * (radius_cm ** 3)
    units["volume_cm3"] = "cm^3"

    mass_g = volume_cm3 * PT_DENSITY_G_PER_CM3
    units["mass_g"] = "g"

    surface_area_m2 = surface_area_nm2 * NM2_TO_M2
    units["surface_area_m2"] = "m^2"

    total_surface_area_m2 += surface_area_m2
    total_mass_g += mass_g

if total_mass_g <= 0:
    raise ValueError("Computed total Pt mass is zero or negative, so surface area per gram cannot be calculated.")

# Apply correction factor for cuboctahedron shape
corrected_total_surface_area_m2 = total_surface_area_m2 * SURFACE_AREA_CORRECTION_FACTOR
specific_surface_area_m2_per_g = corrected_total_surface_area_m2 / total_mass_g
units["specific_surface_area_m2_per_g"] = "m^2/g"

count_difference = particle_count - bin_counts_total
units["count_difference"] = "count"

# Calculate average and standard deviation of particle sizes
average_diameter_nm = statistics.mean(diameters_nm)
standard_deviation_nm = statistics.stdev(diameters_nm)

# Outputs
with open(OUTPUT_FILE, "w", encoding="utf-8") as out_file:
    out_file.write("Particle Size Distribution of Pt Nanoparticles\n")
    out_file.write(f"Input file: {INPUT_FILE_1}\n")
    out_file.write(f"Diameter range: {MIN_DIAMETER_NM:.2f} to {MAX_DIAMETER_NM:.2f} nm\n")
    out_file.write(f"Bin width: {bin_width_nm:.2f} nm\n")
    out_file.write(f"Valid particles analyzed: {particle_count}\n")
    out_file.write(f"Average diameter: {average_diameter_nm:.2f} nm\n")
    out_file.write(f"Standard deviation: {standard_deviation_nm:.2f} nm\n")
    out_file.write("\n")
    out_file.write("Bin_center_nm\tCount\tFrequency\n")

    for row in distribution_rows:
        bin_center_nm_out = row[0]
        units["bin_center_nm_out"] = "nm"

        count_out = row[1]
        units["count_out"] = "count"

        frequency_out = row[2]
        units["frequency_out"] = "fraction"

        out_file.write(f"{bin_center_nm_out:.2f}\t{count_out}\t{frequency_out:.6f}\n")

    out_file.write("\n")
    out_file.write("Theoretical Pt Surface Area Calculation\n")
    out_file.write(f"Pt density (used): {PT_DENSITY_G_PER_CM3:.4f} g/cm^3\n")
    out_file.write(f"Total Pt surface area: {corrected_total_surface_area_m2:.12e} m^2\n")
    out_file.write(f"Total Pt mass: {total_mass_g:.12e} g\n")
    out_file.write(f"Specific Pt surface area: {specific_surface_area_m2_per_g:.6f} m^2/g\n")

print(f"Particle size distribution saved to: {OUTPUT_FILE}")
print(f"Valid particles analyzed: {particle_count}")
print(f"Specific Pt surface area: {specific_surface_area_m2_per_g:.6f} m^2/g")
print(f"Average diameter: {average_diameter_nm:.2f} nm")
print(f"Standard deviation: {standard_deviation_nm:.2f} nm")

print("=== Variable Summary ===")
for k, v in units.items():
    try:
        if not isinstance(eval(k), (list, dict, set)):
            print(f"  {k} = {eval(k)} [{v}]")
    except:
        pass