drop_lib.py

"""Functions for analyzing droplet shape from molecular dynamics simulations.

This module provides routines for sampling and smoothing the droplet
profile, fitting local lines or global circle to determine the contact angle
and contact radius, as well as visualizing the results.

"""
import matplotlib.pyplot as plt
import numpy as np

from typing import Optional
from numpy.typing import ArrayLike

def compute_com(
    atom_positions: np.ndarray, 
    atom_masses: Optional[np.ndarray] = None
) -> np.ndarray:
    """Compute the center-of-mass (COM) of given atoms.
    
    Parameter
    ---------
    atom_positions : ndarray
        Array of shape (N, 3) containing atomic positions.
    atom_masses : ndarray, optional
        Array of (N,) containing masses. 
        If None, all masses are assumed equal.
    
    Returns
    -------
    com : ndarray 
        Length 4 as [cx, cy, cz, cr]. 
        Returns NaN if no atoms are provided.
    
    Notes
    -----
    The center of mass is calculated as:
    
        COMq = sum( qi * mi )/sum( mi )  where qi = (x, y, z) of atom i
    
    """
    if atom_positions.size == 0:
        return np.ones(4) * np.nan
    
    # Constructing an array of equal masses of one:
    if atom_masses is None:
        atom_masses = np.ones_like(atom_positions[:,0])
    
    total_mass = atom_masses.sum()
    com_xyz = atom_positions.T @ atom_masses / total_mass
    com_r = np.linalg.norm(com_xyz)
    return np.concatenate([ com_xyz, [com_r] ])


def compute_rg(
    atom_positions: np.ndarray, 
    atom_masses: Optional[np.ndarray] = None
) -> float:
    """Compute the radius of gyration of atoms about their COM.

    Parameters
    ----------
    atom_positions : ndarray
        Array of shape (N, 3) containing atomic positions.
    atom_masses : ndarray, optional
        Array of (N,) containing masses. 
        If None, all masses are assumed equal.
                    
    Returns
    -------
    rg : float
        Radius of Gyration 
        Returns NaN if no atoms are provided.
        
    Notes
    ----
    The radius of gyration is calculated as:
    
        rg = sqrt( sum(mi * (ri - rcm)**2) / sum( mi ) )
        
    where 'ri' and 'mi' are the position vector and mass of particle i.
    'rcm' is the position vector of COM.
    
    """
    if atom_positions.size == 0:
        return np.nan
    
    # Constructing an array of equal masses of one:
    if atom_masses is None:
        atom_masses = np.ones_like(atom_positions[:,0])
        
    com_xyz = compute_com(atom_positions, atom_masses)[:3]      
    disp = atom_positions - com_xyz   
    rg2 = np.sum(atom_masses * np.sum(disp**2, axis=1)) / atom_masses.sum()
    return float(np.sqrt(rg2))


def fit_taubin(coords: np.ndarray
) -> tuple[float, float, float, float]:
    """Fit a circle to (x, y) points by Taubin's algebraic SVD method.

    Nikolai Chernov's implementation based on Taubin's method.

    References
    ----------
    1. N. Chernov, Circular and Linear Regression
        (doi: 10.1201/EBK1439835906)
    2. G. Taubin, Estimation Of Planar Curves, Surfaces And 
        Nonplanar Space Curves Defined By Implicit Equations.
        (doi: 10.1109/34.103273)
            
    Parameters
    ----------
    coords : ndarray
        (x, y) coordinates.

    Returns
    -------
    xc : float
        x-coordinate of circle center.
    yc : float 
        y-coordinate of circle center.
    r : float
        radius of circle.
    s : float
        RMS error of the fit.

    """
    # ensure input is array
    arr = np.asarray(coords)

    # compute centroid and center data
    centroid = arr.mean(axis=0)
    Xc = arr[:, 0] - centroid[0]
    Yc = arr[:, 1] - centroid[1]

    # compute standardized radial distances
    radial_sq = Xc**2 + Yc**2
    mean_radial = radial_sq.mean()
    norm_radial = (radial_sq - mean_radial) / (2 * np.sqrt(mean_radial))

    # build data matrix and perform SVD
    data_matrix = np.column_stack((norm_radial, Xc, Yc))
    Vt = np.linalg.svd(data_matrix, full_matrices=False)[2]
    V = Vt.T

    # extract and scale coefficients
    coeff = V[:, 2]
    coeff[0] /= 2 * np.sqrt(mean_radial)
    coeff = np.hstack((coeff, -mean_radial * coeff[0]))

    # compute circle center and radius
    offset = -(coeff[1:3]) / (2 * coeff[0])
    xc, yc = offset + centroid
    r = np.sqrt(coeff[1]**2 + coeff[2]**2 - 4 * coeff[0] * coeff[3]) / (2 * abs(coeff[0]))

    # compute RMS error (sigma)
    dx = arr[:, 0] - xc
    dy = arr[:, 1] - yc
    s = np.sqrt(np.mean((np.sqrt(dx**2 + dy**2) - r)**2))

    return xc, yc, r, s


def smooth_mod_sinc(
    data: ArrayLike, 
    ms1: bool = False,
    deg: int = 2,
    half_width: int = 10
) -> np.ndarray:
    """Smooths input data using a modified sinc kernel smoother.

    References
    ----------
    Schmid et al., Why and How Savitzky-Golay Filters Should 
        Be Replaced. (doi: 10.1021/acsmeasuresciau.1c00054)

    Parameters
    ----------
    data : array_like
        1-D sequence of input values to smooth.
    ms1 : bool, optional
        If True enables MS1 variant (default=False).
    deg : int
        Even degree between 2 and 10 (default=2).
    half_width : int
        Kernel half-width (default=10).

    Returns
    -------
    smoothed : ndarray
        1-D array of smoothed `data` .
        
    Raises
    ------
    ValueError
        If `deg` is not an even integer in [2,10], or if `half_width`
        is smaller than the required minimum for the chosen variant.

    """
    # Correction tables for MS0 and MS1 variants
    CORRECTION_MS0 = [
        None, None, None,
        [[0.001717576, 0.02437382, 1.64375]],
        [[0.0043993373, 0.088211164, 2.359375],
         [0.006146815, 0.024715371, 3.6359375]],
        [[0.0011840032, 0.04219344, 2.746875],
         [0.0036718843, 0.12780383, 2.7703125]]
    ]
    CORRECTION_MS1 = [
        None, None,
        [[0.021944195, 0.050284006, 0.765625]],
        [[0.0018977303, 0.008476806, 1.2625],
         [0.023064667, 0.13047926, 1.2265625]],
        [[0.0065903002, 0.057929456, 1.915625],
         [0.0023234477, 0.010298849, 2.2726562],
         [0.021046653, 0.16646601, 1.98125]],
        [[9.749618e-4, 0.0020742896, 3.74375],
         [0.008975366, 0.09902466, 2.7078125],
         [0.0024195414, 0.010064855, 3.296875],
         [0.019185117, 0.18953617, 2.784961]]
    ]

    # Validate parameters
    if deg < 2 or deg > 10 or deg % 2 != 0:
        raise ValueError(f"Invalid deg {deg}; must be even and between 2 and 10.")
    min_hw = (deg // 2 + 1) if ms1 else (deg // 2 + 2)
    if half_width < min_hw:
        raise ValueError(f"Invalid half_width {half_width}; must be >= {min_hw}.")

    # Select correction data and compute coefficients
    corr_table = CORRECTION_MS1 if ms1 else CORRECTION_MS0
    params = corr_table[deg // 2]
    coeffs = None
    if params:
        coeffs = [a + b / (c - half_width) ** 3 for a, b, c in params]

    # Build normalization kernel
    size = half_width + 1
    kernel = np.zeros(size)
    total = 0
    decay = 2 if ms1 else 4
    factor = deg + (2 if ms1 else 4)

    for i in range(size):
        x = i / size
        arg = np.pi * 0.5 * factor * x
        window = 1.0 if i == 0 else np.sin(arg) / arg
        # add correction terms
        if coeffs:
            for idx, coef in enumerate(coeffs):
                window += coef * x * np.sin((2 * idx + 1) * np.pi * x)
        # apply exponential decay
        window *= np.exp(-x * x * decay)

        kernel[i] = window
        total += window if i == 0 else 2 * window

    kernel /= total

    # Prepare boundary extension
    zero_pt = size / (1 + 0.5 * deg) if ms1 else size / (1.5 + 0.5 * deg)
    beta = (0.65 + 0.35 * np.exp(-0.55 * (deg - 4))) if ms1 else (0.70 + 0.14 * np.exp(-0.60 * (deg - 4)))
    L = int(np.ceil(zero_pt * beta))
    win = np.cos(0.5 * np.pi / (zero_pt * beta) * np.arange(L)) ** 2
    
    # Prepare the data
    arr = np.asarray(data)
    n = len(arr)
    ext = np.zeros(n + 2 * half_width)
    ext[half_width:half_width + n] = arr

    # Estimate edge slopes for extension
    L0 = min(L, n)
    s0 = np.polyfit(np.arange(L0), arr[:L0], 1, w=win[:L0])[0]
    s1 = np.polyfit(np.arange(L0), arr[-L0:], 1, w=win[:L0])[0]
    ext[:half_width] = arr[0] - s0 * np.arange(half_width, 0, -1)
    ext[-half_width:] = arr[-1] + s1 * np.arange(1, half_width + 1)

    # Convolve and extract result
    result = np.zeros_like(ext)
    for i in range(half_width, len(ext) - half_width):
        val = kernel[0] * ext[i]
        for j in range(1, size):
            val += kernel[j] * (ext[i - j] + ext[i + j])
        result[i] = val

    return result[half_width:half_width + n]


def extract_profile_2d(
    positions: np.ndarray,
    h_intf: float, 
    include_layer: bool = False,
    separate_sides: bool = False,
    dimensionality: str = '3d'
) -> np.ndarray | tuple[np.ndarray, np.ndarray]:
    """Extract a 2D profile of a sessile droplet.
    
    After recentering on its center of mass, projects a 3D (spherical
    cap) or 2D (cylindrical) droplet into the xz-plane. Molecules up 
    to 1 or 2 factors of `h_intf` are removed.
    
    Parameters
    ----------
    positions : ndarray
        Array of shape (N, 3) of (x, y, z) coordinates.
    h_intf : float
        Height threshold for excluding the interfacial layer along z.
    include_layer : bool, optional 
        If True, excludes points within twice the `h_intf`. 
        False ignores points below interfacial layer.
        (default=False)
    separate_sides : bool, optional
        If True, return two profiles (left, right); otherwise, return
        a single combined profile. (default=False)
    dimensionality : str
        '3d' or '2d'.
        
    Returns
    -------
    profile2d : ndarray or tuple of ndarray
        If `separate_sides` is False, an array of shape (M, 2);
        if True, a tuple `(profile_left, profile_right)`.
            
    Raises:
        ValueError: if `dimensionality` is not '2d' or '3d'.
 
    Notes
    -----
    Radial distances for 3D are computed as:
    
        x_new = sign(x) * sqrt( x**2 + y**2 )
        
    """
    
    if dimensionality not in ['2d', '3d']:
        raise ValueError(f"Unsupported drop_type: {dimensionality!r}."
                         "Use '2d' or '3d'."
                         )
        
    if dimensionality == '3d':  
        # convert to 2D profile
        profile2d = np.zeros_like(positions)
        radial_distances = np.linalg.norm(positions[:, :2], axis=1)
        profile2d[:, 0] = np.where(positions[:, 0] < 0, -radial_distances,
                                   radial_distances)
        profile2d[:, 2] = positions[:, 2]  # Preserve z-coordinates
        
    elif dimensionality == '2d':
        # profile is already 2D:
        profile2d = positions.copy()
    
    # Remove y-coordinates
    profile2d = np.delete(profile2d, 1, axis=1)
    
    # separate left and right profiles
    profile_left = profile2d[profile2d[:, 0] < 0]
    # ensure left x-values are positive (visualization purposes)
    profile_left[:, 0] = np.abs(profile_left[:, 0])
    profile_right = profile2d[profile2d[:, 0] > 0]
    
    # exclude upto 1 or 2 times the interfacial layer's height.
    exclude_height = h_intf if include_layer else 2*h_intf
    profile_right = profile_right[profile_right[:, 1] > exclude_height]
    profile_left = profile_left[profile_left[:, 1] > exclude_height]

    # if no profile atom were extracted, return None
    if (profile_right.size == 0) or (profile_left.size == 0):
        return None, None 

    # return the entire left and right profiles as a single profile:
    if not separate_sides:
        profile_left [:, 0] = -profile_left[:, 0]
        profile2d = np.r_[profile_left, profile_right]
        return profile2d 

    return profile_left, profile_right


def resample_polar(
    x_data: ArrayLike, 
    y_data: ArrayLike, 
    num_points: int = 10
) -> tuple[np.ndarray, np.ndarray]:
    """Resample droplet profile in polar coordinates.
    
    Converts droplet profile data to polar coordinates for improved
    sampling. Sorts profile data according to the polar angle `theta`.
    Evenly samples the `theta` values. Then, interpolates `r` at 
    new `theta` values.
    
    Parameters
    ----------
    x_data : array_like
        1-D input values for x.
    y_data : array_like: 
        1-D input values for y.
    num_points : int
        Number of points to sample along the arc. (default=10)

    Returns
    -------
    r_new : ndarray
        Resampled radii at new theta points (zeros are filtered out).
    theta_new : ndarray
        Corresponding polar angles to `r_new`.

    Raises
    ------
    ValueError
        If `x_data` and `y_data` differ in lenght or have fewer 
        than 2 points.
    
    """
    
    x_data = np.array(x_data, dtype=float)
    y_data = np.array(y_data, dtype=float)
    
    if len(x_data) != len(y_data) or len(x_data) < 2:
        raise ValueError("x_data and y_data must be non-empty"
                         "and of the same length.")
    
    # Convert to polar coordinates
    theta = np.arctan2(y_data, x_data)
    r = np.sqrt(x_data**2 + y_data**2)
    
    # Detect outliers
    msk_out = detect_outliers(r)
    theta = theta[~msk_out]
    r = r[~msk_out]
    
    # Sort by theta
    sort_idx = np.argsort(theta)
    theta_sorted = theta[sort_idx]
    r_sorted = r[sort_idx]
    
    # Identify the min and max angles
    theta_min, theta_max = theta_sorted[0], theta_sorted[-1]
    
    # Define evenly for theta
    theta_new = np.linspace(theta_min, theta_max, num_points)

    # Interpolate r at new theta_nonlinear points
    r_new = np.interp(theta_new, theta_sorted, r_sorted)
    
    # Do not return empty bin samples
    return r_new[r_new!=0], theta_new[r_new!=0]


def detect_outliers(y_data: np.ndarray, z_thresh: int = 2.2) -> np.ndarray:
    """Detects outliers in the polar profile due to thermal fluctuations.
    
    Parameters
    ----------
    y_data : array_like
        1-D array of data for outlier detection.
    z_thresh : int
        Z-score threshold for outliers (default=2.2).
    
    Returns
    -------
    outlier_mask : array_like
        Boolean mask of `y_data`, True for outliers.
    
    Raises
    ------
    ValueError
        If `y_data` is not 1-D.

    Notes
    -----
    Data is split into non-overlapping windows of size N//10. The final
    window, if smaller than `window_size`, is merged with the previous window. 
    Z-scores are computed within each window, and points exceeding `z_thresh` 
    are flagged.
    
    """
    
    arr = np.asarray(y_data, dtype=float)
    if arr.ndim != 1:
        raise ValueError("y_data must be a 1-D array")

    n = arr.shape[0]
    if n < 2:
        return np.zeros(0, dtype=bool)

    window_size = max(1, n // 10)
    outlier_mask = np.zeros(n, dtype=bool)

    start = 0
    while start < n:
        end = start + window_size
        if end > n:
            # Merge final short chunk with the preceding window
            start = max(0, n - window_size)
            end = n

        chunk = arr[start:end]
        mean = chunk.mean()
        std = chunk.std()
        if std and not np.isnan(std):
            zscores = (chunk - mean) / std
            outlier_mask[start:end] = np.abs(zscores) > z_thresh

        if end == n:
            break
        start += window_size

    return outlier_mask


def smooth_profile_polar(
    profile2d: np.ndarray,
    rad_curvature: float,
    num_points: int = 50
) -> tuple[np.ndarray, np.ndarray]:
    """Samples and smooths shifted droplet profile in polar coordinates.
    
    Parameters
    ----------
    profile2d : ndarray
        Shape (N, 2) of (x, y, z) coordinates of the 2D profile.
    rad_curvature : float
        Radius of curvature from a fitted circle.
    num_points : int, optional
        Number of points to sample for unifrom and smoothed profiles.
        (default=50)
    
    Returns
    -------
    profile_uniform: ndarray
        Uniformly sampled droplet profile in Carthesian coordinates
        with shape (num_points, 2)
    profile_smoothed: ndarray
        Smoothed droplet profile in Cartesian coordinates (same shape).
    
    Notes
    -----
    This functions shifts the profile by `rad_curvature` and call
    ``resample_polar`` to evenly sample the radial distances. It then
    ``smooth_mod_sinc`` (MS1=True, half_width=num_points//2) to smooth
    sampled radii before converting back to Cartesian coordinates.
    
    """

    # Get positions
    x_positions, y_positions = profile2d[:, 0], profile2d[:, 1]
    
    # Move y coordinates and resample the proilfe evenly
    r_sample, theta_sample = resample_polar(y_positions - rad_curvature,
                                            x_positions,
                                            num_points)
    
    # Move sampled y coordinates back to original
    y_uniform = r_sample * np.cos(theta_sample) + rad_curvature
    x_uniform = r_sample * np.sin(theta_sample)
    profile_uniform = np.c_[x_uniform, y_uniform]
    
    # Smooth polar radius of the coordinates
    # `deg = 2` and `m = total number of coordinate poitns / 2` recommended
    r_smoothed = smooth_mod_sinc(r_sample,
                                 deg=2,
                                 ms1=True,
                                 half_width=max(1, num_points // 2) )
    
    # Convert smoothed profile to Carthesian coordinates
    # and moves y values back to original
    y_smoothed = r_smoothed * np.cos(theta_sample) + rad_curvature
    x_smoothed = r_smoothed * np.sin(theta_sample)
    profile_smoothed = np.c_[x_smoothed, y_smoothed]
    return profile_uniform, profile_smoothed


def compute_cir_params(
    profile2d: np.ndarray,
    h_intf: float,
    num_points: int = 50
) -> tuple[np.ndarray, np.ndarray, np.ndarray, float, float]:
    """Computes wetting parameters of a droplet using circle fitting.
    
    Parameters
    ----------
    profile2d : np.ndarray
        Array of shape (N, 2) with columns (x, z) of the 2-D profile.
    h_intf : float
        Height threshold for excluding the interfacial layer along z.
    num_points : int, optional
        Number of points to return for the fitted profile(default=50).
        
    Returns
    -------
    polyfit_xy : ndarray      
        (M1, 2) points used for linear fit (for graphing).
    tangent_xy : ndarray      
        (M2, 2) tangent line points (for graphing).
    circle_profile : np.ndarray
        (M3, 2) fitted circle points above `h_intf`.
        Uses `fit_taubin`.
    contact_radius : float
        contact radius extrapolated to the interfacial layer.
    contact_angle : float
        contact angle in degrees by half-angle method.
        
    Notes
    -----
    Fits a circle to (x, z) via `fit_taubin`. Only keeps the arc 
    where `z >= h_intf`. Contact radius is evaluated to `z = h_intf`. 
    Contact angle is computed using the half-angle method
    
        CA = 2 * arctan( CH / CR )
        
    where CH and CR are contact height and contact radius. If the fitted
    circle does not intersect the surface, returns `contact_radius=0` and
    and `contact_angle=180`.

    """
    
    # fit the circle to the profile:
    fit_params = fit_taubin(profile2d)
    thetas = np.linspace(-np.pi, np.pi, num_points*10)

    # circle points:
    x_circ = fit_params[0] + np.cos(thetas)*fit_params[2]
    y_circ = fit_params[1] + np.sin(thetas)*fit_params[2]
    circle_profile = np.c_[x_circ[y_circ >= h_intf], 
                           y_circ[y_circ >= h_intf]]
    
    # determine contact with the substrate:
    dy =  h_intf - fit_params[1] 
    if np.abs(dy) > fit_params[2]:     
        return (np.zeros((1, 2)), np.zeros((1, 2)), circle_profile, 0., 180.)

    contact_height = y_circ.max() - h_intf
    dx = np.sqrt(fit_params[2]**2 - dy**2)
    contact_radius = dx + fit_params[0]

    # contact angle by half-angle method
    contact_angle = np.degrees( 2 * np.arctan2(contact_height, contact_radius) )

    # the rest are purely for graphing purposes:
    # mask data up to twice interfacial height.   
    msk4graph = np.where( (y_circ <= 2*h_intf)
                        & (y_circ >= h_intf)
                        & (x_circ > 0))[0]
    x4poly = x_circ[msk4graph]
    y4poly = y_circ[msk4graph]
    polyfit_xy = np.c_[x4poly, y4poly] 
    prm = np.polyfit(y4poly, x4poly, 1)

    # Points to graph for tangent line:
    y = np.arange(h_intf, x4poly.max()/2)
    tangent_xy = np.c_[np.polyval(prm, y), y]

    return (polyfit_xy, tangent_xy, circle_profile, contact_radius, 
            contact_angle)


def compute_local_params(
    smoothed_profile: np.ndarray,
    h_intf: float,
    z_first_peak: Optional[float] = None,
    z_second_peak: Optional[float] = None
) -> tuple[np.ndarray, np.ndarray, float, float]:
    """Computes wetting parameters of a droplet using local fitting.

    Parameters
    ----------
    smoothed_profile : ndarray
        Array of shape (N, 2) with columns (x, y) of smoothed 2-D profile.
    h_intf : float
        Height threshold for excluding the interfacial layer along z.
    z_first_peak : float 
        Height of the first density peak. If None, defaults to `0.6*h_intf`.
    z_second_peak : float
        Height of the second density peak. If None, defaults to `2*h_intf`. 

    Returns
    -------
    polyfit_xy : ndarray      
        (M1, 2) points used for local linear fit (for graphing).
    tangent_xy : ndarray      
        (M2, 2) tangent line points (for graphing).
    contact_radius : float
        contact radius extrapolated to `z = z_first_peak`.
    contact_angle : float
        contact angle in degrees by local-fit method.
        
    Notes
    -----
    Extracts point below `z_second_peak` in the `smoothed_profile`.
    Fits a line to the points and reports use its slope to compute the
    `contact_angle`. `contact_radius` is determined by extrapolating the
    tangent to `z = z_first_peak`.
    
    Raises
    ------
    ValueError
        The `z_first_peak` must be smaller than `z_second_peak`.
    """
    
    # Given `h_intf`, 0.6 is a safe choice for the peak position
    if z_first_peak is None:
        z_first_peak = 0.6*h_intf
    
    # Given `h_intf`, 2 is a safe choice for the peak position
    if z_second_peak is None:
        z_second_peak = 2*h_intf
        
    if z_second_peak <= z_first_peak:
        raise ValueError(" 'z_first_peak' must be < 'z_second_peak' ")
        
    # Extract data and fit tangent line
    x4poly = smoothed_profile[:, 0]
    y4poly = smoothed_profile[:, 1]
    
    # Mask data up to twice interfacial width to correct for CL bending
    msk_2nd_layer = np.where(y4poly < z_second_peak)[0]
    x4poly = smoothed_profile[msk_2nd_layer, 0]
    y4poly = smoothed_profile[msk_2nd_layer, 1]
    
    # in case of early droplet contact with the substrate, 
    # no points would be detected:
    if y4poly.shape[0] < 2:
        return (np.zeros((0, 2)), np.zeros((0, 2)), 0., 180.)
    
    polyfit_xy = np.c_[x4poly, y4poly]  
    prm = np.polyfit(y4poly, x4poly, 1)
    
    # Contact angle based on the slope of the line
    p_slope = np.polyder(prm, 1)
    contact_angle = 90 + float(np.degrees(np.arctan2(p_slope[0], 1.0)))

    # Points to graph tangent line:
    y = np.arange(z_first_peak, x4poly.max()/2)
    tangent_xy = np.c_[np.polyval(prm, y), y]    
    
    # Contact radius as the first point in the fitted line
    if tangent_xy.shape[0] == 0:
        contact_radius = 0.
    else:
        contact_radius = tangent_xy[0, 0]
            
    return polyfit_xy, tangent_xy, contact_radius, contact_angle


def analyze_droplet(
    intf_molecs: np.ndarray, 
    z_interfacial_layer: float, 
    z_surface: float = 0,
    z_peak1: Optional[float] = None, 
    z_peak2: Optional[float] = None,
    side: str = 'left',
    dim: str = '3d',
    visualize: bool = False
) -> tuple[float, float]:
    """Analyzes the droplet profile to characterize its wetting behavior.
    
    Parameters
    ----------
    intf_molecs : ndarray
        Array of shape (N, 3) with (x, y, z) interfacial molecule positions.
    z_interfacial_layer : float
        Height threshold for excluding the interfacial layer along z.
    z_surface : float
        Location of substrate along z (default=0).
    z_peak1 : float, optional
        Height of the first density peak. If None, `compute_local_params`
        defaults to `0.6*h_intf`.
    z_peak2 : float, optional
        Height of the second density peak. If None, `compute_local_params`
        defaults to `2*h_intf`.
    side : {'left', 'right', 'both'}, optional
        Determines which side of the droplet is analyzed by 
        `compute_local_params`. If `both`, defaults to circle fitting via
        `compute_cir_params` (default=`left`).
    dim : {'2d', '3d'}, optional
        Droplet dimensionality. Cylindrical (`2d`) and spherical (`3d`)
        (default=`3d`).
    visualize : str, optional
        If True, visualizes the profiles and fits (default=False).
    
    Returns
    -------
    contact_angle : float
        Contact angle in degress.
    contact_radius : float
        Contact radius in x units.
        
    Notes
    -----
    Translates the interfacial points based on the radius of curvature of the droplet. 
    Transforms the positions in cylindrical coordinates to extract a two-dimensional
    profile using `extract_profile_2d`. Invokes `compute_local_params` for `side=left`
    and `side=right` and uses `compute_cir_params` for `side=both`.
    
    Raises
    ------
    ValueError
        If `side` is not one of {`left`, `right`, `both`} or if `dim` is not
        one of {`2d`, `3d`}
    """
    
    if side not in ['left', 'right', 'both']:
        raise ValueError(f"Unsupported analysis: {side!r}"
                         "Use only `left`, `right`, or `both`.")
    
    if dim not in ['2d', '3d']:
        raise ValueError(f"Unsupported dimensions: {dim!r}"
                         "Use  `2d` for cylindrical and "
                         "`3d` for spherical droplets.")
    
    intf_molecs = intf_molecs.copy()
    # Recenter the droplet
    intf_molecs[:, :2] -= compute_com(intf_molecs)[:2] 
    # Place surface at zero
    intf_molecs[:, 2] -= z_surface
    
    # First, the profile is shifted by `y_correction` (profile's radius of curvature) for better sampling:
    drop_profile =  extract_profile_2d(intf_molecs, z_interfacial_layer, include_layer=False, separate_sides=False, dimensionality=dim)
    
    if drop_profile[0] is None:
        return 0. , 180.
    
    # Fit a circle to the droplet for correction value
    y_correction = fit_taubin(drop_profile)[1]
    
    # Determine the optimal number of points for smoothing and fitting
    opt_num_pts = int( 0.5 * drop_profile[:, :2].max())    # Flexibility analysis should be run on this value.

    # Extract left and right profiles:
    profile_left, profile_right =  extract_profile_2d(
        intf_molecs, z_interfacial_layer, 
        include_layer=True, separate_sides=True, dimensionality=dim)
    
    if isinstance(drop_profile, tuple) and drop_profile[0] is None:
        return 0.0, 180.0
    
    if (profile_left is None) or (profile_right is None):
        return 0. , 180.
    
    if side.lower() == 'left':
        # Grab left-side of the profile
        drop_profile = profile_left.copy()
        # Smooth the profile
        uniform_profile, smoothed_profile =  smooth_profile_polar(drop_profile, y_correction, num_points=opt_num_pts)
        # Finally, we calculate the local parameters
        poly_xy, tang_xy, contact_radius, contact_angle = compute_local_params(smoothed_profile, z_interfacial_layer, z_first_peak=z_peak1, z_second_peak=z_peak2)       
       
    elif side.lower() == 'right':
        # Grab right-side of the profile
        drop_profile = profile_right.copy()
        # Smooth the profile
        uniform_profile, smoothed_profile =  smooth_profile_polar(drop_profile, y_correction, num_points=opt_num_pts)
        # Finally, we calculate the local parameters
        poly_xy, tang_xy, contact_radius, contact_angle = compute_local_params(smoothed_profile, z_interfacial_layer, z_first_peak=z_peak1, z_second_peak=z_peak2)
        
    elif side.lower() == 'both':
        drop_profile = extract_profile_2d(intf_molecs, z_interfacial_layer, include_layer=True, separate_sides=False, dimensionality=dim)
        poly_xy, tang_xy, uniform_profile, contact_radius, contact_angle = compute_cir_params(drop_profile, z_interfacial_layer, num_points=opt_num_pts)
        smoothed_profile = uniform_profile
    
    if (contact_radius < 0) or (drop_profile is None):
        return 0., 180.

    # Simple visualization to check the fitting procedure:
    if visualize:
        graph_droplet_analysis(drop_profile, uniform_profile, smoothed_profile, tang_xy, poly_xy, contact_radius, z_interfacial_layer, analysis_side=side)
    
    return float(contact_radius), float(contact_angle)


def graph_droplet_analysis(
    drop_profile: np.ndarray,
    uniform_profile: np.ndarray, 
    smoothed_profile: np.ndarray,
    tangent_line: np.ndarray,
    points4tangent: np.ndarray,
    contact_radius: float,
    z_interfacial_layer: float,
    analysis_side: str='left'
) -> None:
    """Visualizes the analysis, extracted profiles, and computed fits.
    
    Parameters
    ----------
    drop_profile : ndarray
        (N, 2) array of original profile points.
    uniform_profile : ndarray
        (M1, 2) array of uniformly sampled profile points.
    smoothed_profile : ndarray
        (M2, 2) array of smoothed profile points. If `analysis_side=both`
        this is not plotted.
    tangent_line : ndarray
        (M3, 2) points of the tangent line visualizing the contact angle.
    points4tangent : ndarray
        (M4, 2) array of points used for contact angle calculation in 
        `compute_local_params`. If `analysis_side=both`, purely for
        visualization purposes.
    contact_radius : float
        Contact radius value.
    z_interfacial_layer : float
        Height threshold for excluding the interfacial layer along z.
    analysis_side : str
        If one of {`left`, `right`}, plots the `uniform_profile`.
        If `both`, skips plotting (default=`left`)
        
    Returns
    -------
    None
    
    """
    
    _, ax = plt.subplots(figsize=(4, 4))
    # Original droplet profile
    ax.plot(drop_profile[:, 0], drop_profile[:, 1],
             'C0.', label="Original", alpha=0.1)
    
    # Smoothed or fitted profile
    if analysis_side == 'both':
        ax.plot(smoothed_profile[:, 0], smoothed_profile[:, 1],
                'C3.', label="Fitted", lw=2, zorder=99)
    else:
        ax.plot(smoothed_profile[:, 0], smoothed_profile[:, 1],
                'C3.', label="Smoothed", lw=2, zorder=99)
    # contact radius
    ax.plot([0, contact_radius],
            [z_interfacial_layer/2, z_interfacial_layer/2], 'kd--', ms=5)
    # points used for fitting the tangent
    ax.plot(points4tangent[:, 0], points4tangent[:, 1],
             'C1o', label="Sampled")
    ax.set_xlabel(r"$x-\rm{axis~(\AA)}$", fontsize=14)
    ax.set_ylabel(r"$z-\rm{axis~(\AA)}$", fontsize=14)
    # slope at three-phase point
    ax.plot(tangent_line[:, 0], tangent_line[:, 1],
             'k-', lw=2, label='Slope', zorder=100)
    
    if analysis_side in ['left', 'right']:
        # Uniformly sampled profile
        ax.plot(uniform_profile[:, 0], uniform_profile[:, 1],
                 'C0', lw=2, label="Uniform")
    
    max_dim = np.max([drop_profile[:, 0].max(), drop_profile[:, 1].max()])
    ax.set_xlim(0, max_dim + 10)
    ax.set_ylim(0, max_dim + 10)
    ax.legend(loc='upper right')
    ax.set_title(f'Analysis for the {analysis_side} side')
    plt.show()
    plt.close()