data.py

#!/usr/bin/env python3
"""
Summary plots
"""

import sys
import math
import pandas as pd
import matplotlib.pyplot as plt


STEADY_STATE_FRACTION = 0.30
PUBLICATION_MODE = True


def apply_publication_style():
    if not PUBLICATION_MODE:
        return
    plt.rcParams.update({
        "font.size": 14,
        "axes.titlesize": 15,
        "axes.labelsize": 14,
        "xtick.labelsize": 12,
        "ytick.labelsize": 12,
        "legend.fontsize": 12,
        "figure.titlesize": 15,
        "lines.linewidth": 1.8,
    })


def ivantsov_2d(P):
    try:
        import mpmath as mp
    except Exception:
        return None

    if P <= 0.0:
        return 0.0
    s = mp.sqrt(P)
    return float(mp.sqrt(mp.pi * P) * mp.e**P * mp.erfc(s))


def clean_series(s):
    s = pd.to_numeric(s, errors="coerce")
    s = s.replace([math.inf, -math.inf], math.nan)
    return s


def trim_startup(df, vel_col=None):
    out = df.copy()
    if vel_col is not None and vel_col in out.columns:
        v = clean_series(out[vel_col])
        vmax = v.max(skipna=True)
        if pd.notna(vmax) and vmax > 0.0:
            keep = v >= 0.02 * vmax
            if keep.any():
                first = keep.idxmax()
                return out.loc[first:].reset_index(drop=True)
    return out


def col_lookup(df, name):
    lname = name.lower()
    for c in df.columns:
        if c.lower() == lname:
            return c
    return None


def col_any(df, *names):
    for n in names:
        c = col_lookup(df, n)
        if c is not None:
            return c
    return None

def print_series_summary(label, s):
    vals = clean_series(s)
    vals = vals[vals.notna()]
    if len(vals) == 0:
        return
    mean = vals.mean()
    std = vals.std()
    print(f"{label:<20s}: {mean:.5f} ± {std:.5f}")


def new_figure():
    plt.figure(figsize=(7.2, 5.2))


def plot_ca_series(ax, x, y, label="CA", color="tab:blue", marker="o"):
    """Helper to consistently plot CA data with visible points only."""
    m = x.notna() & y.notna()
    if not m.any():
        return

    ax.plot(
        x[m], y[m],
        marker=marker,
        linestyle="None",
        markersize=5.5,
        markeredgecolor="k",
        markeredgewidth=0.6,
        label=label,
        color=color
    )


def main():
    if len(sys.argv) < 2:
        print("Usage: python data3.py data.csv")
        return 2

    apply_publication_style()

    path = sys.argv[1]
    df = pd.read_csv(path)
    df.columns = [c.strip() for c in df.columns]

    c_time  = col_lookup(df, "time(s)")
    c_vel   = col_lookup(df, "vel(um/s)")
    c_ucool = col_lookup(df, "undercoolavg(degC)")
    c_rad   = col_lookup(df, "radiusKGT(um)")
    c_sol   = col_any(df, "solcap", "DFRZ")
    c_p     = col_lookup(df, "peclet")
    c_om    = (col_lookup(df, "Omega") or col_lookup(df, "OmegaCA"))
    c_omkgt = (col_lookup(df, "OmegaKGT") or
               col_lookup(df, "OmegaKGT(2D)") or
               col_lookup(df, "OmegaKGT(3D)"))

    c_dtkgt = col_lookup(df, "DTKGT(degC)")
    c_cl    = col_lookup(df, "Clstar(wt%)")
    c_theta = col_lookup(df, "thetaBOX(deg)")

    sol_label = "SOLCAP" if col_lookup(df, "solcap") else "DFRZ"

    required = [c_time, c_vel, c_ucool, c_rad, c_sol, c_p, c_om]
    if any(c is None for c in required):
        print("Missing one or more required lean CSV columns.")
        return 1

    df = df.sort_values(c_time).reset_index(drop=True)
    df_trim = trim_startup(df, c_vel)

    # 1. Omega vs Peclet (with clear CA points)
    x = clean_series(df_trim[c_p])
    y = clean_series(df_trim[c_om])
    if x.notna().any() and y.notna().any():
        new_figure()
        ax = plt.gca()
        plot_ca_series(ax, x, y, label="CA data", color="tab:blue")

        # Ivantsov reference
        try:
            pmin = float(max(1.0e-6, x[x > 0].min()))
            pmax = float(max(pmin * 1.001, x.max()))
            lpmin = math.log(pmin)
            lpmax = math.log(pmax)
            curve_p = []
            curve_o = []
            for i in range(200):
                a = i / 199.0
                pp = math.exp(lpmin + a * (lpmax - lpmin))
                vv = ivantsov_2d(pp)
                if vv is not None:
                    curve_p.append(pp)
                    curve_o.append(vv)
            if curve_p:
                plt.plot(curve_p, curve_o, label="2-D Ivantsov", color="tab:orange", linewidth=2.2)
        except Exception:
            pass

        plt.xlabel("Peclet number, P = V R / (2 D)")
        plt.ylabel("Supersaturation, Omega")
        plt.title("Omega vs Peclet")
        plt.legend()
        plt.grid(True)
        plt.tight_layout()

    # 2. Supersaturation vs time
    x = clean_series(df_trim[c_time])
    y = clean_series(df_trim[c_om])
    if x.notna().any() and y.notna().any():
        new_figure()
        ax = plt.gca()
        plot_ca_series(ax, x, y, label="CA", color="tab:blue")

        if c_omkgt is not None:
            y2 = clean_series(df_trim[c_omkgt])
            m2 = x.notna() & y2.notna()
            if m2.any():
                plt.plot(x[m2], y2[m2], label="KGT", color="tab:orange", linewidth=2.0)

        plt.xlabel("Time (s)")
        plt.ylabel("Supersaturation, Omega")
        plt.title("Supersaturation vs time")
        plt.legend()
        plt.grid(True)
        plt.tight_layout()

    # 3. Velocity vs time
    x = clean_series(df_trim[c_time])
    y = clean_series(df_trim[c_vel])
    if x.notna().any() and y.notna().any():
        new_figure()
        ax = plt.gca()
        plot_ca_series(ax, x, y, label="CA velocity", color="tab:blue")
        plt.xlabel("Time (s)")
        plt.ylabel("Velocity (um/s)")
        plt.title("Tip velocity vs time")
        plt.grid(True)
        plt.tight_layout()

    # 4. Undercooling vs time
    x = clean_series(df_trim[c_time])
    y = clean_series(df_trim[c_ucool])
    if x.notna().any() and y.notna().any():
        new_figure()
        ax = plt.gca()
        plot_ca_series(ax, x, y, label="CA", color="tab:blue")

        if c_dtkgt is not None:
            y2 = clean_series(df_trim[c_dtkgt])
            m2 = x.notna() & y2.notna()
            if m2.any():
                plt.plot(x[m2], y2[m2], label="KGT", color="tab:orange", linewidth=2.0)

        plt.xlabel("Time (s)")
        plt.ylabel("Undercooling (degC)")
        plt.title("Tip undercooling vs time")
        plt.legend()
        plt.grid(True)
        plt.tight_layout()

    # 5. KGT radius vs time
    x = clean_series(df_trim[c_time])
    y = clean_series(df_trim[c_rad])
    if x.notna().any() and y.notna().any():
        new_figure()
        ax = plt.gca()
        plot_ca_series(ax, x, y, label="CA", color="tab:blue")
        plt.xlabel("Time (s)")
        plt.ylabel("KGT tip radius (um)")
        plt.title("KGT tip radius vs time")
        plt.grid(True)
        plt.tight_layout()

    # 6. SOLCAP vs time
    x = clean_series(df_trim[c_time])
    y = clean_series(df_trim[c_sol])
    if x.notna().any() and y.notna().any():
        new_figure()
        ax = plt.gca()
        plot_ca_series(ax, x, y, label="CA", color="tab:blue")
        plt.xlabel("Time (s)")
        plt.ylabel(sol_label)
        plt.title(f"{sol_label} vs time")
        plt.grid(True)
        plt.tight_layout()

    # Optional: Interface composition
    if c_cl is not None:
        x = clean_series(df_trim[c_time])
        y = clean_series(df_trim[c_cl])
        if x.notna().any() and y.notna().any():
            new_figure()
            ax = plt.gca()
            plot_ca_series(ax, x, y, label="CA", color="tab:blue")
            plt.xlabel("Time (s)")
            plt.ylabel("Clstar (wt%)")
            plt.title("Interface liquid composition vs time")
            plt.grid(True)
            plt.tight_layout()

    # Optional: BOX orientation
    if c_theta is not None:
        x = clean_series(df_trim[c_time])
        y = clean_series(df_trim[c_theta])
        if x.notna().any() and y.notna().any():
            new_figure()
            ax = plt.gca()
            plot_ca_series(ax, x, y, label="CA", color="tab:blue")
            plt.xlabel("Time (s)")
            plt.ylabel("thetaBOX (deg)")
            plt.title("BOX orientation vs time")
            plt.grid(True)
            plt.tight_layout()

    # Steady-state summary
    n = len(df_trim)
    if n > 10:
        i0 = int((1.0 - STEADY_STATE_FRACTION) * n)
        df_ss = df_trim.iloc[i0:].copy()

        print()
        print("===============================================")
        print("LEAN CSV STEADY-STATE SUMMARY")
        print("===============================================")
        print(f"Points used           : {len(df_ss)}")
        print_series_summary("Velocity (um/s)", df_ss[c_vel])
        print_series_summary("Undercooling (degC)", df_ss[c_ucool])
        print_series_summary("RadiusKGT (um)", df_ss[c_rad])
        print_series_summary("SOLCAP", df_ss[c_sol])
        print_series_summary("Peclet", df_ss[c_p])
        print_series_summary("Omega", df_ss[c_om])
        if c_omkgt is not None:
            print_series_summary("OmegaKGT", df_ss[c_omkgt])
        if c_dtkgt is not None:
            print_series_summary("DTKGT (degC)", df_ss[c_dtkgt])
        if c_cl is not None:
            print_series_summary("Clstar (wt%)", df_ss[c_cl])
        if c_theta is not None:
            print_series_summary("thetaBOX (deg)", df_ss[c_theta])
        print("===============================================")

    plt.show()
    return 0


if __name__ == "__main__":
    raise SystemExit(main())