mcltspice/src/mcp_ltspice/power_analysis.py

"""Power and efficiency calculations from simulation data."""

import numpy as np

# np.trapz was renamed to np.trapezoid in numpy 2.0
_trapz = getattr(np, "trapezoid", getattr(np, "trapz", None))


def compute_instantaneous_power(voltage: np.ndarray, current: np.ndarray) -> np.ndarray:
    """Element-wise power: P(t) = V(t) * I(t).

    Args:
        voltage: Voltage waveform array
        current: Current waveform array

    Returns:
        Instantaneous power array (same length as inputs)
    """
    return np.real(voltage) * np.real(current)


def compute_average_power(time: np.ndarray, voltage: np.ndarray, current: np.ndarray) -> float:
    """Time-averaged power: P_avg = (1/T) * integral(V*I dt).

    Uses trapezoidal integration over the full time span.

    Args:
        time: Time array in seconds
        voltage: Voltage waveform array
        current: Current waveform array

    Returns:
        Average power in watts
    """
    t = np.real(time)
    duration = t[-1] - t[0]

    if len(t) < 2 or duration <= 0:
        return 0.0

    p_inst = np.real(voltage) * np.real(current)
    return float(_trapz(p_inst, t) / duration)


def compute_efficiency(
    time: np.ndarray,
    input_voltage: np.ndarray,
    input_current: np.ndarray,
    output_voltage: np.ndarray,
    output_current: np.ndarray,
) -> dict:
    """Compute power conversion efficiency.

    Args:
        time: Time array in seconds
        input_voltage: Input voltage waveform
        input_current: Input current waveform
        output_voltage: Output voltage waveform
        output_current: Output current waveform

    Returns:
        Dict with efficiency_percent, input_power_watts,
        output_power_watts, power_dissipated_watts
    """
    p_in = compute_average_power(time, input_voltage, input_current)
    p_out = compute_average_power(time, output_voltage, output_current)
    p_dissipated = p_in - p_out

    if abs(p_in) < 1e-15:
        efficiency = 0.0
    else:
        efficiency = (p_out / p_in) * 100.0

    return {
        "efficiency_percent": efficiency,
        "input_power_watts": p_in,
        "output_power_watts": p_out,
        "power_dissipated_watts": p_dissipated,
    }


def compute_power_spectrum(
    time: np.ndarray,
    voltage: np.ndarray,
    current: np.ndarray,
    max_harmonics: int = 20,
) -> dict:
    """FFT of instantaneous power to identify ripple frequencies.

    Args:
        time: Time array in seconds
        voltage: Voltage waveform array
        current: Current waveform array
        max_harmonics: Maximum number of frequency bins to return

    Returns:
        Dict with frequencies, power_magnitudes, dominant_freq, dc_power
    """
    t = np.real(time)

    if len(t) < 2:
        return {
            "frequencies": [],
            "power_magnitudes": [],
            "dominant_freq": 0.0,
            "dc_power": 0.0,
        }

    dt = (t[-1] - t[0]) / (len(t) - 1)
    if dt <= 0:
        return {
            "frequencies": [],
            "power_magnitudes": [],
            "dominant_freq": 0.0,
            "dc_power": float(np.mean(np.real(voltage) * np.real(current))),
        }

    p_inst = np.real(voltage) * np.real(current)
    n = len(p_inst)

    spectrum = np.fft.rfft(p_inst)
    freqs = np.fft.rfftfreq(n, d=dt)
    magnitudes = np.abs(spectrum) * 2.0 / n
    magnitudes[0] /= 2.0  # DC component correction

    dc_power = float(magnitudes[0])

    # Dominant AC frequency (largest non-DC bin)
    if len(magnitudes) > 1:
        dominant_idx = int(np.argmax(magnitudes[1:])) + 1
        dominant_freq = float(freqs[dominant_idx])
    else:
        dominant_freq = 0.0

    # Trim to requested harmonics (plus DC)
    limit = min(max_harmonics + 1, len(freqs))
    freqs = freqs[:limit]
    magnitudes = magnitudes[:limit]

    return {
        "frequencies": freqs.tolist(),
        "power_magnitudes": magnitudes.tolist(),
        "dominant_freq": dominant_freq,
        "dc_power": dc_power,
    }


def compute_power_metrics(time: np.ndarray, voltage: np.ndarray, current: np.ndarray) -> dict:
    """Comprehensive power report for a voltage/current pair.

    Power factor here is defined as the ratio of real (average) power
    to apparent power (Vrms * Irms).

    Args:
        time: Time array in seconds
        voltage: Voltage waveform array
        current: Current waveform array

    Returns:
        Dict with avg_power, rms_power, peak_power, min_power, power_factor
    """
    v = np.real(voltage)
    i = np.real(current)
    p_inst = v * i

    if len(p_inst) == 0:
        return {
            "avg_power": 0.0,
            "rms_power": 0.0,
            "peak_power": 0.0,
            "min_power": 0.0,
            "power_factor": 0.0,
        }

    avg_power = compute_average_power(time, voltage, current)
    rms_power = float(np.sqrt(np.mean(p_inst**2)))
    peak_power = float(np.max(p_inst))
    min_power = float(np.min(p_inst))

    # Power factor = avg power / apparent power (Vrms * Irms)
    v_rms = float(np.sqrt(np.mean(v**2)))
    i_rms = float(np.sqrt(np.mean(i**2)))
    apparent = v_rms * i_rms

    if apparent < 1e-15:
        power_factor = 0.0
    else:
        power_factor = avg_power / apparent

    return {
        "avg_power": avg_power,
        "rms_power": rms_power,
        "peak_power": peak_power,
        "min_power": min_power,
        "power_factor": power_factor,
    }