dna/chem5/plot.py

0001 import numpy as np
0002 import matplotlib.pyplot as plt
0003 import pandas as pd
0004
0005 # Load data from Species.txt
0006 fileName1 = "Species.txt"
0007 data = pd.read_csv(fileName1, delim_whitespace=True, header=None, names=["Time", "Value", "Err", "Species"],
0008                    comment="#")
0009
0010 # Load data from e_aq.txt
0011 fileName2 = "e_aq.txt"
0012 data_eaq = pd.read_csv(fileName2, delim_whitespace=True, header=None, names=["Time", "Value", "Err"], comment="#")
0013 data_eaq["Species"] = "e_aq^-1"  # Assign species name
0014
0015 # Load data from OH.txt
0016 fileName3 = "OH.txt"
0017 data_oh = pd.read_csv(fileName3, delim_whitespace=True, header=None, names=["Time", "Value", "Err"], comment="#")
0018 data_oh["Species"] = "°OH^0"  # Assign species name
0019
0020 # Log scale transformation (avoid log(0) error)
0021 data["Time"] = np.log10(data["Time"].replace(0, np.nan))
0022 data_eaq["Time"] = np.log10(data_eaq["Time"].replace(0, np.nan))
0023 data_oh["Time"] = np.log10(data_oh["Time"].replace(0, np.nan))
0024
0025 # Define species to plot with LaTeX formatting
0026 species_list = ["°OH^0", "e_aq^-1", "H3O^1", "H2O2^0", "H_2^0", "H^0"]
0027 labels = [
0028     r"$\mathrm{\cdot OH}$",  # OH radical
0029     r"$\mathrm{e^-_{aq}}$",  # Aqueous electron
0030     r"$\mathrm{H_3O^+}$",  # Hydronium ion
0031     r"$\mathrm{H_2O_2}$",  # Hydrogen peroxide
0032     r"$\mathrm{H_2}$",  # Molecular hydrogen
0033     r"$\mathrm{H^\bullet}$"  # Hydrogen radical
0034 ]
0035 marker_styles = ["s", "^"]  # Square (e_aq), Triangle (OH)
0036
0037 # Create figure and axes
0038 fig, axes = plt.subplots(2, 3, figsize=(14, 8), sharex=True)
0039 axes = axes.flatten()
0040
0041 for i, (species, label) in enumerate(zip(species_list, labels)):
0042     ax = axes[i]
0043     subset = data[data["Species"] == species]
0044
0045     if not subset.empty:
0046         ax.plot(subset["Time"], subset["Value"], linestyle='-', marker='o', label="Simulation", color='black')
0047         ax.errorbar(subset["Time"], subset["Value"], yerr=subset["Err"], fmt='o', capsize=5, color='black')
0048
0049     # Add e_aq.txt as additional points
0050     if species == "e_aq^-1":
0051         ax.scatter(data_eaq["Time"], data_eaq["Value"], marker=marker_styles[0], color='red',
0052                    label=r"$\mathrm{e^-_{aq} \ exp}$")
0053
0054     # Add OH.txt as additional points
0055     if species == "°OH^0":
0056         ax.scatter(data_oh["Time"], data_oh["Value"], marker=marker_styles[1], color='blue',
0057                    label=r"$\mathrm{\cdot OH \ exp}$")
0058
0059     ax.set_title(label, fontsize=14)
0060     ax.set_yscale("linear")
0061     ax.grid(True)
0062     ax.legend()
0063
0064 # Set common labels
0065 fig.text(0.5, 0.001, r"$\mathrm{Time \ (log(ps))}$", ha='center', fontsize=14)
0066 fig.text(0.005, 0.5, r"$\mathrm{G(Species/100 \ eV)}$", va='center', rotation='vertical', fontsize=14)
0067
0068 plt.tight_layout()
0069 plt.show()