COMSOL

#!/usr/bin/env python3
import sys
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from io import StringIO
from scipy.interpolate import griddata

def parse_comsol_complex(s):
    """
    Convert a COMSOL-style complex string (e.g. '1.2E-5-3.4E-6i')
    into a Python complex number.
    Returns 0j if s is blank or '0'.
    """
    s = s.strip()
    if s == '' or s == '0':
        return 0 + 0j
    s = s.replace('i', 'j')
    return complex(s)

def main(doPlot = True, csv_file = R"C:\Users\winga\Desktop\UCLA\Research\COMSOL\15Pert\37.955GHz O MODE\COMSOL_OUTPUT_37.9GHZ_15.0000.csv",freq = "37.955",x1 = 2.3644,x2=2.404,y1=-0.070778,y2=0.070778):

    # -------------------------------
    # Step 1: Read and parse the COMSOL CSV file
    # -------------------------------
    with open(csv_file, "r", encoding="utf-8") as f:
        lines = f.readlines()

    metadata_lines = []
    header_line = None
    data_lines = []
    for line in lines:
        stripped = line.strip()
        if not stripped:
            continue
        if header_line is None:
            if stripped.startswith('% X'):
                header_line = stripped
            else:
                metadata_lines.append(stripped)
        else:
            data_lines.append(stripped)

    # (Optional) Parse metadata
    expected_keys = ["Model", "Version", "Date", "Dimension",
                     "Nodes", "Expressions", "Description", "Length unit"]
    metadata = {}
    for meta_line in metadata_lines:
        if meta_line.startswith('%'):
            clean_line = meta_line.lstrip('%').strip()
            parts = clean_line.split(',', 1)
            if len(parts) == 2:
                key = parts[0].strip()
                value = parts[1].strip()
                if key in expected_keys:
                    metadata[key] = value
    print("Parsed Metadata:")
    for key in expected_keys:
        print(f"{key}: {metadata.get(key, 'Not Found')}")

    if header_line is None:
        raise ValueError("No data header line found (starting with '% X').")
    clean_header = header_line.lstrip('%').strip()
    header_cols = clean_header.split(',')
    print("\nData columns (from header):", header_cols)

    combined_csv = StringIO()
    combined_csv.write(",".join(header_cols) + "\n")
    for dl in data_lines:
        combined_csv.write(dl + "\n")
    combined_csv.seek(0)

    df = pd.read_csv(combined_csv, sep=",")
    print("\nDataFrame head:")
    print(df.head())

    # -------------------------------
    # Step 2: Compute the E-field magnitude
    # -------------------------------
    ex_col = "emw.Ex (V/m) @ freq="+freq
    ey_col = "emw.Ey (V/m) @ freq="+freq
    ez_col = "emw.Ez (V/m) @ freq="+freq
    df["Ex"] = df[ex_col].apply(parse_comsol_complex)
    df["Ey"] = df[ey_col].apply(parse_comsol_complex)
    df["Ez"] = df[ez_col].apply(parse_comsol_complex)

    df["E_mag"] = np.sqrt(np.abs(df["Ex"])**2 +
                          np.abs(df["Ey"])**2 )
                          #np.abs(df["Ez"])**2)
    print("E_mag: min =", np.min(df["E_mag"]), "max =", np.max(df["E_mag"]))

    # -------------------------------
    # Step 3: Interpolate scattered E-field data onto a regular 2D grid
    # -------------------------------
    points = df[['X', 'Y']].values
    values = df['E_mag'].values

    x_min, x_max = df['X'].min(), df['X'].max()
    y_min, y_max = df['Y'].min(), df['Y'].max()
    print(f"X domain: {x_min} to {x_max} m")
    print(f"Y domain: {y_min} to {y_max} m")

    # Use a moderate grid resolution for memory efficiency:
    n_grid_x = 3000
    n_grid_y = 3000
    print(f"Using grid resolution: {n_grid_x} x {n_grid_y}")

    grid_x = np.linspace(x_min, x_max, n_grid_x)
    grid_y = np.linspace(y_min, y_max, n_grid_y)
    grid_X, grid_Y = np.meshgrid(grid_x, grid_y)

    E2D = griddata(points, values, (grid_X, grid_Y), method='linear')
    E2D = np.nan_to_num(E2D)
    print("Interpolated grid shape:", E2D.shape)
    print("Interpolated E2D: min =", np.min(E2D), "max =", np.max(E2D))

    # -------------------------------
    # Step 4: Define the port and an offset line for integration
    # -------------------------------
    # Given port endpoints:
    xp1 = x1    # 2.3423 m
    yp1 = y1 # -0.0466 m
    xp2 = x2    # 2.3645 m
    yp2 = y2  # 0.0466 m
    p1 = np.array([xp1, yp1])
    p2 = np.array([xp2, yp2])

    port_vec = p2 - p1
    L_port = np.linalg.norm(port_vec)
    print(f"Port length = {L_port:.4f} m")

    # Unit vector along the port:
    u_port = port_vec / L_port
    # Choose a normal that points in the direction of the backscattered beam.
    # Previously we set: u_norm = [u_port[1], -u_port[0]]
    u_norm = np.array([u_port[1], -u_port[0]])
    # If needed, flip sign: u_norm = -u_norm

    # Define an offset distance d_offset (in m) to shift the integration line behind the port.
    d_offset = 0.005  # e.g. 5 mm behind the port
    p1_offset = p1 + d_offset * u_norm
    p2_offset = p2 + d_offset * u_norm

    # Parameterize the offset line:
    npts_line = 200
    s_array = np.linspace(0, 1, npts_line)
    line_points = np.array([p1_offset + s*(p2_offset - p1_offset) for s in s_array])

    # -------------------------------
    # Step 5: Perform a line integral along the offset line
    # -------------------------------
    # For a free-space wave, the time-averaged Poynting flux is:
    # S = 0.5 * eps0 * c * |E|^2   (W/m²)
    # In 2D, if the simulation is per unit out-of-plane length, then integrating S along a line gives the total power (W).
    eps0 = 8.854187817e-12  # F/m
    c = 3.0e8               # m/s

    # For each point on the line, interpolate E2D from the grid:
    E_line = griddata(points, values, (line_points[:,0], line_points[:,1]), method='linear')
    E_line = np.nan_to_num(E_line)
    S_line = 0.5 * eps0 * c * (E_line**2)

    # The differential element along the line:
    dl = L_port / (npts_line - 1)
    total_power_line = np.sum(S_line) * dl
    print(f"\nIntegrated power crossing the offset line: {total_power_line:.3e} W")

    # -------------------------------
    # Step 6: Plot the results
    # -------------------------------
    # Plot the interpolated E-field heatmap and overlay the port and offset line.
    width = x_max - x_min
    height = y_max - y_min
    aspect_ratio = width / height

    plt.figure(figsize=(8,6))
    im = plt.imshow(E2D, origin="lower",
                    extent=[x_min, x_max, y_min, y_max],
                    cmap="jet", vmin=np.min(E2D), vmax=np.max(E2D))
    plt.gca().set_aspect(aspect_ratio)
    plt.colorbar(im, label="|E| (V/m)")
    plt.xlabel("X [m]")
    plt.ylabel("Y [m]")
    plt.title("Backscattered Electric Field Magnitude")

    # Draw the original port line and the offset line
    plt.plot([p1[0], p2[0]], [p1[1], p2[1]], 'w-', linewidth=3, label="Port")
    plt.plot([p1_offset[0], p2_offset[0]], [p1_offset[1], p2_offset[1]],
             'w--', linewidth=3, label="Offset integration line")

    plt.legend()
    plt.savefig('lineint.jpg',dpi=300)
    if doPlot:
        plt.show()

    # Additionally, plot the Poynting flux along the offset line
    plt.figure(figsize=(8,5))
    plt.plot(s_array * L_port, S_line, 'bo-', linewidth=2, markersize=4)
    plt.xlabel("Distance along offset line [m]")
    plt.ylabel("Poynting flux S [W/m²]")
    plt.title("Poynting Flux along the Offset Integration Line")
    plt.grid(True)
    plt.savefig('pynting.jpg',dpi=300)

    if doPlot:
        plt.show()
    return total_power_line

if __name__ == "__main__":
    main()