Atomizer/docs/plans/DYNAMIC_RESPONSE_IMPLEMENTATION_PLAN.md

# Dynamic Response Processor - Implementation Plan

## Atomizer Integration Assessment & Revised Roadmap

**Document Version**: 1.1
**Created**: 2025-12-22
**Revised**: 2025-12-22
**Status**: Ready for Implementation
**Based On**: DYNAMIC_RESPONSE_MASTER_PLAN.md

---

## Executive Summary

After analyzing the Master Plan against Atomizer's current architecture, this document provides a **refined implementation plan** that:

1. Aligns with established patterns (extractors, insights, registry)
2. Leverages existing code (modal_mass extractor, modal_analysis insight)
3. Creates the new `processors/` layer for algorithm code
4. Integrates seamlessly with dashboard and protocols

**Key Finding**: The Master Plan is well-designed. This document reorganizes it to fit Atomizer's conventions and identifies reuse opportunities.

---

## 1. Architecture Alignment

### Current Atomizer Layers

```
┌─────────────────────────────────────────────────────────────┐
│                    USER INTERFACE                            │
│  (Dashboard / CLI / run_optimization.py)                     │
└─────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌─────────────────────────────────────────────────────────────┐
│                    EXTRACTORS                                │
│  optimization_engine/extractors/extract_*.py                 │
│  - Simple function interface: extract_X(file) → dict         │
│  - Called from objective functions                           │
│  - Example: extract_frequency(), extract_modal_mass()        │
└─────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌─────────────────────────────────────────────────────────────┐
│                    PROCESSORS (NEW)                          │
│  optimization_engine/processors/dynamic_response/            │
│  - Algorithm implementations (modal superposition, FRF)      │
│  - Class-based with rich state                               │
│  - Called by extractors                                      │
└─────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌─────────────────────────────────────────────────────────────┐
│                    INSIGHTS                                  │
│  optimization_engine/insights/*.py                           │
│  - Visualization classes with @register_insight              │
│  - Generate HTML/Plotly for dashboard                        │
│  - Example: ModalInsight, StressFieldInsight                 │
└─────────────────────────────────────────────────────────────┘
```

### Proposed Integration

```
optimization_engine/
├── extractors/
│   ├── __init__.py                    # Add new exports
│   ├── extract_modal_mass.py          # EXISTING - F06 parsing ✓
│   ├── extract_random_response.py     # NEW - High-level random vib interface
│   ├── extract_sine_response.py       # NEW - Sine sweep interface
│   └── extract_dynamic_stress.py      # NEW - Dynamic stress interface
│
├── processors/                         # NEW DIRECTORY
│   ├── __init__.py
│   └── dynamic_response/
│       ├── __init__.py                # Public API
│       ├── modal_database.py          # ModeShape, ModalDatabase
│       ├── transfer_functions.py      # TransferFunctionEngine
│       ├── random_vibration.py        # RandomVibrationProcessor
│       ├── sine_sweep.py              # SineSweepProcessor (Phase 4)
│       ├── psd_profiles.py            # PSD definitions
│       └── utils/
│           ├── __init__.py
│           ├── frequency.py           # Frequency array utilities
│           ├── integration.py         # Spectral integration
│           └── statistics.py          # Peak factors
│
├── insights/
│   ├── __init__.py                    # Add new registration
│   ├── modal_analysis.py              # EXISTING - Mode shapes ✓
│   └── dynamic_response.py            # NEW - Grms, PSD plots
│
└── hooks/
    └── dynamic_response/               # NEW
        ├── pre_analysis.py            # Validate modal data
        └── post_analysis.py           # Cache results
```

---

## 2. Code Reuse Opportunities

### Existing Assets to Leverage

| Asset | Location | Reuse Strategy |
|-------|----------|----------------|
| Modal mass extraction | `extract_modal_mass.py` | Call from `ModalDatabase.from_op2_f06()` |
| F06 MEFFMASS parsing | `_parse_modal_effective_mass()` | Import directly |
| Participation factors | `_parse_participation_factors()` | Import directly |
| OP2 reading | `pyNastran` (installed) | Use for eigenvectors |
| Modal insight base | `modal_analysis.py` | Extend or reference patterns |
| Insight registry | `@register_insight` decorator | Use for DynamicResponseInsight |
| Extractor patterns | All `extract_*.py` | Follow return dict pattern |

### Key Code Connections

```python
# In processors/dynamic_response/modal_database.py
from optimization_engine.extractors.extract_modal_mass import (
    extract_modal_mass,
    _parse_modal_effective_mass,  # Internal but useful
)

class ModalDatabase:
    @classmethod
    def from_op2_f06(cls, op2_file, f06_file, **options):
        # Use existing extractor for participation factors
        modal_data = extract_modal_mass(f06_file, mode=None)

        # Extract eigenvectors from OP2
        from pyNastran.op2.op2 import OP2
        op2 = OP2()
        op2.read_op2(str(op2_file))

        # Build ModeShape objects combining both sources
        ...
```

---

## 3. Implementation Phases (Revised)

### Phase 1: Core Infrastructure (3-4 days)

**Goal**: Create processor foundation with test-validated components.

#### 1.1 Create Directory Structure
```bash
mkdir -p optimization_engine/processors/dynamic_response/utils
touch optimization_engine/processors/__init__.py
touch optimization_engine/processors/dynamic_response/__init__.py
touch optimization_engine/processors/dynamic_response/utils/__init__.py
```

#### 1.2 Implement Data Classes

**File**: `optimization_engine/processors/dynamic_response/modal_database.py`

```python
"""
Modal Database - Central storage for modal analysis data.

Combines:
- Eigenvectors from OP2 (pyNastran)
- Participation factors from F06 (existing extract_modal_mass)
- User-specified damping

Provides efficient access for dynamic calculations.
"""

from dataclasses import dataclass, field
from typing import Dict, List, Optional, Tuple
from pathlib import Path
import numpy as np
import logging

logger = logging.getLogger(__name__)


@dataclass
class ModeShape:
    """Single mode shape container."""

    mode_number: int
    frequency: float  # Hz
    damping_ratio: float = 0.02

    # Modal properties (from F06)
    generalized_mass: float = 1.0
    participation_x: float = 0.0
    participation_y: float = 0.0
    participation_z: float = 0.0
    effective_mass_x: float = 0.0
    effective_mass_y: float = 0.0
    effective_mass_z: float = 0.0

    # Eigenvector at nodes (from OP2)
    # Dict[node_id, np.array([tx, ty, tz, rx, ry, rz])]
    eigenvector: Dict[int, np.ndarray] = field(default_factory=dict)

    # Modal stress shapes (optional, for stress recovery)
    # Dict[element_id, np.array([sxx, syy, szz, sxy, syz, sxz, vm])]
    stress_shapes: Dict[int, np.ndarray] = field(default_factory=dict)

    @property
    def omega(self) -> float:
        """Angular frequency (rad/s)."""
        return 2 * np.pi * self.frequency

    @property
    def omega_squared(self) -> float:
        return self.omega ** 2

    @property
    def quality_factor(self) -> float:
        """Q = 1/(2ζ)"""
        return 1 / (2 * self.damping_ratio) if self.damping_ratio > 0 else np.inf

    def get_participation(self, direction: str) -> float:
        """Get participation factor for direction ('x', 'y', 'z')."""
        return getattr(self, f'participation_{direction}', 0.0)

    def get_effective_mass(self, direction: str) -> float:
        return getattr(self, f'effective_mass_{direction}', 0.0)


class ModalDatabase:
    """
    Central modal data storage with lazy loading.

    Combines OP2 (eigenvectors) and F06 (participation) data.
    Provides caching via HDF5 for optimization loops.
    """

    def __init__(
        self,
        modes: List[ModeShape],
        total_mass: float,
        model_name: str = "unnamed"
    ):
        self.modes = modes
        self.total_mass = total_mass
        self.model_name = model_name
        self._cache = {}

    @property
    def n_modes(self) -> int:
        return len(self.modes)

    @property
    def frequency_range(self) -> Tuple[float, float]:
        if not self.modes:
            return (0.0, 0.0)
        freqs = [m.frequency for m in self.modes]
        return (min(freqs), max(freqs))

    def get_frequencies(self) -> np.ndarray:
        return np.array([m.frequency for m in self.modes])

    def get_participation_factors(self, direction: str) -> np.ndarray:
        return np.array([m.get_participation(direction) for m in self.modes])

    def get_effective_masses(self, direction: str) -> np.ndarray:
        return np.array([m.get_effective_mass(direction) for m in self.modes])

    def get_cumulative_mass_fraction(self, direction: str) -> np.ndarray:
        """Cumulative effective mass fraction for mode selection."""
        eff_masses = self.get_effective_masses(direction)
        return np.cumsum(eff_masses) / self.total_mass

    def modes_for_mass_fraction(self, fraction: float, direction: str) -> int:
        """Number of modes needed to capture given mass fraction."""
        cumulative = self.get_cumulative_mass_fraction(direction)
        indices = np.where(cumulative >= fraction)[0]
        return int(indices[0] + 1) if len(indices) > 0 else self.n_modes

    @classmethod
    def from_op2_f06(
        cls,
        op2_file: Path,
        f06_file: Optional[Path] = None,
        damping: float = 0.02,
        node_ids: Optional[List[int]] = None,
        element_ids: Optional[List[int]] = None,
    ) -> 'ModalDatabase':
        """
        Build ModalDatabase from NX Nastran results.

        Args:
            op2_file: Path to OP2 with eigenvectors
            f06_file: Path to F06 with MEFFMASS (auto-detected if None)
            damping: Default damping ratio for all modes
            node_ids: Only extract eigenvectors at these nodes
            element_ids: Only extract stress shapes for these elements
        """
        from pyNastran.op2.op2 import OP2
        from optimization_engine.extractors.extract_modal_mass import extract_modal_mass

        op2_path = Path(op2_file)

        # Auto-detect F06
        if f06_file is None:
            f06_path = op2_path.with_suffix('.f06')
            if not f06_path.exists():
                f06_path = None
        else:
            f06_path = Path(f06_file)

        # Extract participation factors from F06
        modal_mass_data = None
        if f06_path and f06_path.exists():
            result = extract_modal_mass(f06_path, mode=None)
            if result['success']:
                modal_mass_data = result.get('modes', [])

        # Read OP2 for eigenvectors
        logger.info(f"Reading OP2: {op2_path}")
        op2 = OP2()
        op2.read_op2(str(op2_path))

        if not op2.eigenvectors:
            raise ValueError(f"No eigenvectors found in {op2_path}")

        # Get total mass
        total_mass = 1.0  # Default
        if hasattr(op2, 'grid_point_weight') and op2.grid_point_weight:
            # Try to get from GPWG
            for gpwg in op2.grid_point_weight.values():
                total_mass = float(gpwg.mass[0])
                break

        # Build ModeShape objects
        modes = []

        for key, eig in op2.eigenvectors.items():
            n_modes = len(eig.modes) if hasattr(eig, 'modes') else 0

            for i in range(n_modes):
                mode_num = int(eig.modes[i])
                freq = float(eig.cycles[i]) if hasattr(eig, 'cycles') else 0.0

                # Get eigenvector data
                eigenvector = {}
                if node_ids is not None:
                    for nid in node_ids:
                        if nid in eig.node_gridtype:
                            idx = np.where(eig.node_gridtype[:, 0] == nid)[0]
                            if len(idx) > 0:
                                eigenvector[nid] = eig.data[i, idx[0], :]
                else:
                    # Store all nodes (memory intensive for large models)
                    for j, (nid, _) in enumerate(eig.node_gridtype):
                        eigenvector[int(nid)] = eig.data[i, j, :]

                # Merge with F06 data
                participation = {'x': 0.0, 'y': 0.0, 'z': 0.0}
                effective = {'x': 0.0, 'y': 0.0, 'z': 0.0}
                gen_mass = 1.0

                if modal_mass_data and mode_num <= len(modal_mass_data):
                    mm = modal_mass_data[mode_num - 1]
                    participation = {
                        'x': mm.get('participation_x', 0.0) or 0.0,
                        'y': mm.get('participation_y', 0.0) or 0.0,
                        'z': mm.get('participation_z', 0.0) or 0.0,
                    }
                    effective = {
                        'x': mm.get('mass_x', 0.0) or 0.0,
                        'y': mm.get('mass_y', 0.0) or 0.0,
                        'z': mm.get('mass_z', 0.0) or 0.0,
                    }

                mode_shape = ModeShape(
                    mode_number=mode_num,
                    frequency=freq,
                    damping_ratio=damping,
                    generalized_mass=gen_mass,
                    participation_x=participation['x'],
                    participation_y=participation['y'],
                    participation_z=participation['z'],
                    effective_mass_x=effective['x'],
                    effective_mass_y=effective['y'],
                    effective_mass_z=effective['z'],
                    eigenvector=eigenvector,
                )
                modes.append(mode_shape)

            break  # Only process first subcase

        logger.info(f"Built ModalDatabase: {len(modes)} modes, "
                   f"freq range {modes[0].frequency:.1f}-{modes[-1].frequency:.1f} Hz")

        return cls(modes=modes, total_mass=total_mass, model_name=op2_path.stem)

    def summary(self) -> str:
        """Human-readable summary."""
        lines = [
            f"ModalDatabase: {self.model_name}",
            f"  Modes: {self.n_modes}",
            f"  Frequency range: {self.frequency_range[0]:.1f} - {self.frequency_range[1]:.1f} Hz",
            f"  Total mass: {self.total_mass:.3f} kg",
        ]
        return "\n".join(lines)
```

#### 1.3 Implement Transfer Functions

**File**: `optimization_engine/processors/dynamic_response/transfer_functions.py`

```python
"""
Transfer Function Engine - FRF computation from modal data.

Computes:
- Single-mode FRF (SDOF transfer function)
- Multi-mode combination (modal superposition)
- Base excitation transmissibility
"""

import numpy as np
from typing import Optional, Tuple
from .modal_database import ModalDatabase, ModeShape


class TransferFunctionEngine:
    """
    Frequency response function computation engine.

    Uses modal superposition to compute dynamic response
    transfer functions from modal analysis data.
    """

    def __init__(self, modal_db: ModalDatabase):
        self.modal_db = modal_db

    def single_mode_frf(
        self,
        frequencies: np.ndarray,
        mode: ModeShape,
        output_type: str = 'displacement'
    ) -> np.ndarray:
        """
        Compute FRF for a single mode.

        H(f) = 1 / [(ω_n² - ω²) + j·2·ζ·ω_n·ω]

        Args:
            frequencies: Frequency array (Hz)
            mode: ModeShape object
            output_type: 'displacement', 'velocity', or 'acceleration'

        Returns:
            Complex FRF array
        """
        omega = 2 * np.pi * frequencies
        omega_n = mode.omega
        zeta = mode.damping_ratio

        # Base displacement FRF
        H = 1.0 / ((omega_n**2 - omega**2) + 2j * zeta * omega_n * omega)

        if output_type == 'velocity':
            H = H * (1j * omega)
        elif output_type == 'acceleration':
            H = H * (-omega**2)

        return H

    def base_excitation_frf(
        self,
        frequencies: np.ndarray,
        direction: str = 'z',
        output_type: str = 'acceleration',
        node_id: Optional[int] = None,
        n_modes: Optional[int] = None,
    ) -> np.ndarray:
        """
        Compute FRF for base excitation (seismic/vibration).

        Combines modal contributions:
        H_total(f) = Σᵢ Lᵢ · φᵢ · Hᵢ(f) / mᵢ

        Args:
            frequencies: Frequency array (Hz)
            direction: Excitation direction ('x', 'y', 'z')
            output_type: 'displacement', 'velocity', 'acceleration'
            node_id: Specific output node (None = use participation only)
            n_modes: Number of modes (None = 95% mass fraction)

        Returns:
            Complex FRF (transmissibility for acceleration output)
        """
        # Determine number of modes
        if n_modes is None:
            n_modes = self.modal_db.modes_for_mass_fraction(0.95, direction)
        n_modes = min(n_modes, self.modal_db.n_modes)

        omega = 2 * np.pi * frequencies
        H_total = np.zeros(len(frequencies), dtype=complex)

        for i in range(n_modes):
            mode = self.modal_db.modes[i]

            # Participation factor for this direction
            L = mode.get_participation(direction)
            if L == 0:
                continue

            # Mode shape at output location
            if node_id is not None and node_id in mode.eigenvector:
                dir_idx = {'x': 0, 'y': 1, 'z': 2}.get(direction, 2)
                phi = mode.eigenvector[node_id][dir_idx]
            else:
                # Use participation as proxy for overall response
                phi = 1.0

            # Single mode FRF
            H_mode = self.single_mode_frf(frequencies, mode, 'displacement')

            # Add modal contribution
            H_total += L * phi * H_mode / mode.generalized_mass

        # Convert to requested output type
        if output_type == 'velocity':
            H_total = H_total * (1j * omega)
        elif output_type == 'acceleration':
            # For base excitation: transmissibility = 1 + ω²·H_disp
            H_total = 1.0 + (-omega**2) * H_total

        return H_total

    def transmissibility(
        self,
        frequencies: np.ndarray,
        direction: str = 'z',
        n_modes: Optional[int] = None,
    ) -> Tuple[np.ndarray, np.ndarray]:
        """
        Compute acceleration transmissibility magnitude and phase.

        Returns:
            (magnitude, phase_degrees)
        """
        H = self.base_excitation_frf(
            frequencies, direction, 'acceleration', n_modes=n_modes
        )
        return np.abs(H), np.angle(H, deg=True)
```

#### 1.4 Implement PSD Profiles

**File**: `optimization_engine/processors/dynamic_response/psd_profiles.py`

(Use content from Master Plan Section 2 - already well-designed)

#### 1.5 Create Tests

**File**: `tests/test_dynamic_response/test_transfer_functions.py`

```python
"""
Transfer Function Tests

Validates FRF computations against analytical solutions.
"""

import numpy as np
import pytest
from optimization_engine.processors.dynamic_response.modal_database import ModeShape, ModalDatabase
from optimization_engine.processors.dynamic_response.transfer_functions import TransferFunctionEngine


def test_sdof_resonance_amplification():
    """At resonance, |H| ≈ Q for lightly damped SDOF."""
    mode = ModeShape(
        mode_number=1,
        frequency=100.0,  # Hz
        damping_ratio=0.02,  # 2% damping → Q = 25
    )
    db = ModalDatabase(modes=[mode], total_mass=1.0)
    engine = TransferFunctionEngine(db)

    # Evaluate at resonance
    f_res = np.array([100.0])
    H = engine.single_mode_frf(f_res, mode, 'displacement')

    # Expected: |H| = 1/(2*zeta*omega_n^2) at resonance
    omega_n = 2 * np.pi * 100
    expected_mag = 1 / (2 * 0.02 * omega_n**2)

    assert np.isclose(np.abs(H[0]), expected_mag, rtol=0.01)


def test_transmissibility_unity_low_freq():
    """Transmissibility → 1 at low frequencies (rigid body)."""
    mode = ModeShape(mode_number=1, frequency=100.0, damping_ratio=0.02)
    db = ModalDatabase(modes=[mode], total_mass=1.0)
    engine = TransferFunctionEngine(db)

    low_freq = np.array([1.0])  # Well below resonance
    mag, _ = engine.transmissibility(low_freq)

    assert np.isclose(mag[0], 1.0, rtol=0.1)


def test_transmissibility_isolation_high_freq():
    """Transmissibility → 0 at high frequencies (isolation)."""
    mode = ModeShape(mode_number=1, frequency=100.0, damping_ratio=0.02)
    mode.participation_z = 1.0
    mode.effective_mass_z = 1.0
    db = ModalDatabase(modes=[mode], total_mass=1.0)
    engine = TransferFunctionEngine(db)

    high_freq = np.array([1000.0])  # Well above resonance
    mag, _ = engine.transmissibility(high_freq, direction='z')

    # Should be much less than 1
    assert mag[0] < 0.1
```

---

### Phase 2: Random Vibration Processor (3-4 days)

**Goal**: Complete random vibration response computation with Miles' equation validation.

#### 2.1 RandomVibrationProcessor

**File**: `optimization_engine/processors/dynamic_response/random_vibration.py`

```python
"""
Random Vibration Processor

Computes Grms, peak response, and stress from PSD input using modal superposition.
"""

from dataclasses import dataclass
from typing import Dict, Any, Optional, List, Tuple
import numpy as np
from scipy.integrate import trapezoid

from .modal_database import ModalDatabase
from .transfer_functions import TransferFunctionEngine
from .psd_profiles import PSDProfile


@dataclass
class RandomVibrationResult:
    """Container for random vibration results."""
    rms_acceleration: float  # Grms
    peak_acceleration: float  # G (k-sigma)
    rms_displacement: float  # mm
    peak_displacement: float  # mm
    dominant_mode: int  # Mode with highest contribution
    modal_contributions: Dict[int, float]  # Mode → % contribution
    response_psd: Tuple[np.ndarray, np.ndarray]  # (frequencies, psd)


class RandomVibrationProcessor:
    """
    Random vibration response via modal superposition.

    Replaces SOL 111/112 with analytical computation.
    """

    def __init__(
        self,
        modal_db: ModalDatabase,
        default_damping: float = 0.02
    ):
        self.modal_db = modal_db
        self.default_damping = default_damping
        self.frf_engine = TransferFunctionEngine(modal_db)

    def compute_response(
        self,
        psd_profile: PSDProfile,
        direction: str = 'z',
        node_id: Optional[int] = None,
        n_modes: Optional[int] = None,
        frequency_resolution: int = 500,
        peak_sigma: float = 3.0,
        output_type: str = 'acceleration'
    ) -> RandomVibrationResult:
        """
        Compute full random vibration response.

        Args:
            psd_profile: Input PSD specification
            direction: Excitation direction ('x', 'y', 'z')
            node_id: Output location (None = overall)
            n_modes: Number of modes (None = auto 95% mass)
            frequency_resolution: Points in frequency array
            peak_sigma: Sigma level for peak (3.0 = 3-sigma)
            output_type: 'acceleration' or 'displacement'

        Returns:
            RandomVibrationResult with all metrics
        """
        # Generate frequency array (log spacing recommended for PSD)
        f_min = max(psd_profile.freq_min, 1.0)
        f_max = psd_profile.freq_max
        frequencies = np.logspace(
            np.log10(f_min),
            np.log10(f_max),
            frequency_resolution
        )

        # Compute transfer function
        H = self.frf_engine.base_excitation_frf(
            frequencies, direction, output_type, node_id, n_modes
        )

        # Interpolate input PSD to our frequency array
        input_psd = psd_profile.interpolate(frequencies)

        # Response PSD: S_y = |H|² × S_x
        response_psd = np.abs(H)**2 * input_psd

        # RMS via integration
        mean_square = trapezoid(response_psd, frequencies)
        rms = np.sqrt(mean_square)
        peak = peak_sigma * rms

        # Also compute displacement if acceleration was requested
        if output_type == 'acceleration':
            H_disp = self.frf_engine.base_excitation_frf(
                frequencies, direction, 'displacement', node_id, n_modes
            )
            resp_psd_disp = np.abs(H_disp)**2 * input_psd
            rms_disp = np.sqrt(trapezoid(resp_psd_disp, frequencies))
            peak_disp = peak_sigma * rms_disp

            rms_accel = rms
            peak_accel = peak
        else:
            rms_disp = rms
            peak_disp = peak
            rms_accel = 0.0
            peak_accel = 0.0

        # Modal contribution analysis
        contributions = self._compute_modal_contributions(
            frequencies, input_psd, direction, n_modes
        )

        dominant_mode = max(contributions, key=contributions.get) if contributions else 1

        return RandomVibrationResult(
            rms_acceleration=rms_accel,
            peak_acceleration=peak_accel,
            rms_displacement=rms_disp * 1000,  # Convert to mm
            peak_displacement=peak_disp * 1000,
            dominant_mode=dominant_mode,
            modal_contributions=contributions,
            response_psd=(frequencies, response_psd),
        )

    def _compute_modal_contributions(
        self,
        frequencies: np.ndarray,
        input_psd: np.ndarray,
        direction: str,
        n_modes: Optional[int]
    ) -> Dict[int, float]:
        """Compute contribution of each mode to total response."""
        if n_modes is None:
            n_modes = self.modal_db.modes_for_mass_fraction(0.95, direction)
        n_modes = min(n_modes, self.modal_db.n_modes)

        contributions = {}
        total = 0.0

        for i in range(n_modes):
            mode = self.modal_db.modes[i]
            L = mode.get_participation(direction)

            # Single mode response
            H_mode = self.frf_engine.single_mode_frf(frequencies, mode, 'acceleration')
            resp_psd = np.abs(H_mode * L)**2 * input_psd
            mode_ms = trapezoid(resp_psd, frequencies)

            contributions[mode.mode_number] = mode_ms
            total += mode_ms

        # Normalize to percentages
        if total > 0:
            contributions = {k: v/total * 100 for k, v in contributions.items()}

        return contributions

    def compute_grms(
        self,
        psd_profile: PSDProfile,
        direction: str = 'z',
        node_id: Optional[int] = None
    ) -> float:
        """Quick Grms extraction for optimization objective."""
        result = self.compute_response(psd_profile, direction, node_id)
        return result.rms_acceleration

    def compute_peak_g(
        self,
        psd_profile: PSDProfile,
        direction: str = 'z',
        sigma: float = 3.0
    ) -> float:
        """Quick peak G extraction for constraint."""
        result = self.compute_response(psd_profile, direction, peak_sigma=sigma)
        return result.peak_acceleration


def miles_equation(f_n: float, zeta: float, psd_level: float) -> float:
    """
    Miles' equation for SDOF with white noise.

    Grms = sqrt(π/2 × f_n × Q × PSD)

    Args:
        f_n: Natural frequency (Hz)
        zeta: Damping ratio
        psd_level: PSD level (G²/Hz)

    Returns:
        Grms (G)
    """
    Q = 1 / (2 * zeta)
    return np.sqrt(np.pi / 2 * f_n * Q * psd_level)
```

#### 2.2 Tests for Random Vibration

```python
def test_miles_equation_validation():
    """Random vibration Grms matches Miles' equation for SDOF."""
    f_n = 100.0
    zeta = 0.02
    psd_level = 0.04  # G²/Hz

    # Analytical
    grms_miles = miles_equation(f_n, zeta, psd_level)

    # Via processor
    mode = ModeShape(
        mode_number=1,
        frequency=f_n,
        damping_ratio=zeta,
        participation_z=1.0,
        effective_mass_z=1.0,
    )
    db = ModalDatabase(modes=[mode], total_mass=1.0)

    psd = PSDProfile.flat(psd_level, freq_min=10, freq_max=1000)
    processor = RandomVibrationProcessor(db)
    grms_processor = processor.compute_grms(psd, 'z')

    assert np.isclose(grms_processor, grms_miles, rtol=0.05)  # 5% tolerance
```

---

### Phase 3: Atomizer Integration (2-3 days)

**Goal**: Connect to extractors, insights, and dashboard.

#### 3.1 Create High-Level Extractor

**File**: `optimization_engine/extractors/extract_random_response.py`

(Use content from Master Plan with adjustments to import from processors/)

#### 3.2 Update Extractor Registry

**File**: `optimization_engine/extractors/__init__.py` (add to existing)

```python
# Phase 5: Dynamic Response (2025-12-XX)
from .extract_random_response import (
    extract_random_response,
    extract_random_stress,
    get_grms,
    get_peak_g,
    get_peak_stress,
)

# Add to __all__
__all__ += [
    'extract_random_response',
    'extract_random_stress',
    'get_grms',
    'get_peak_g',
    'get_peak_stress',
]
```

#### 3.3 Create Dynamic Response Insight

**File**: `optimization_engine/insights/dynamic_response.py`

```python
"""
Dynamic Response Insight

Visualizes random vibration and sine sweep response data.
"""

from pathlib import Path
from typing import Optional
import numpy as np

from .base import StudyInsight, InsightConfig, InsightResult, register_insight


@register_insight
class DynamicResponseInsight(StudyInsight):
    """
    Dynamic response visualization.

    Shows:
    - Input vs Output PSD (log-log)
    - Transmissibility magnitude/phase
    - Modal contribution bar chart
    - Grms history across trials
    """

    insight_type = "dynamic_response"
    name = "Dynamic Response Analysis"
    description = "Random vibration and frequency response visualization"
    category = "dynamics"
    applicable_to = ["random_vib", "dynamics", "modal", "vibration"]
    required_files = ["*.op2"]

    # ... implementation following existing insight patterns
```

---

### Phase 4: Sine Sweep & Advanced Features (Future)

Defer to after Phase 3 is validated:
- Sine sweep processor
- Shock response spectrum
- Multi-axis excitation
- Fatigue from random

---

## 4. Protocol Documentation

### New Protocol: SYS_17_DYNAMIC_RESPONSE.md

Create `docs/protocols/system/SYS_17_DYNAMIC_RESPONSE.md`:

```markdown
# SYS_17: Dynamic Response Processor

## Overview
Modal superposition-based dynamic response for optimization.

## Capabilities
- Random vibration (Grms, peak acceleration, stress)
- Sine sweep (magnitude, phase, resonance search)
- Shock response spectrum (future)

## Quick Reference

| Function | Description | Output |
|----------|-------------|--------|
| `extract_random_response()` | Full random vib metrics | dict |
| `get_grms()` | Quick Grms for objective | float |
| `get_peak_g()` | Quick peak G for constraint | float |

## Example Usage

```python
from optimization_engine.extractors import extract_random_response

# In optimization objective
result = extract_random_response(
    op2_file='model.op2',
    psd_profile='GEVS',
    direction='z',
    damping=0.02
)

objective_value = result['grms']
constraint_value = result['peak_acceleration']
```

## PSD Profiles Available
- GEVS (NASA standard)
- MIL-STD-810G variants (CAT1, CAT4, CAT7, CAT24)
- ESA PSS-01-401 (qualification, acceptance)
- Custom from CSV

## Validation
- Miles' equation (SDOF) - <1% error
- SOL 111 comparison - <5% error target
```

### Update SYS_12 Extractor Library

Add entries for new extractors:

| ID | Name | Function | Input | Output |
|----|------|----------|-------|--------|
| E23 | Random Response | `extract_random_response()` | .op2/.f06 + PSD | dict |
| E24 | Dynamic Stress | `extract_random_stress()` | .op2/.f06 + PSD | MPa |
| E25 | Quick Grms | `get_grms()` | .op2 + PSD | G |

---

## 5. Implementation Checklist

### Phase 1 Checklist
- [ ] Create `optimization_engine/processors/` directory structure
- [ ] Implement `modal_database.py` with ModeShape + ModalDatabase
- [ ] Implement `transfer_functions.py` with TransferFunctionEngine
- [ ] Implement `psd_profiles.py` with PSDProfile + STANDARD_PROFILES
- [ ] Create `utils/` with frequency, integration, statistics helpers
- [ ] Write unit tests for FRF (SDOF validation)
- [ ] Test ModalDatabase.from_op2_f06() with real NX data

### Phase 2 Checklist
- [ ] Implement `random_vibration.py` with RandomVibrationProcessor
- [ ] Implement modal contribution analysis
- [ ] Write Miles' equation validation test
- [ ] Add stress recovery (optional for Phase 2)
- [ ] Integration test with sample study

### Phase 3 Checklist
- [ ] Create `extract_random_response.py` in extractors
- [ ] Update `extractors/__init__.py` with exports
- [ ] Create `DynamicResponseInsight` in insights
- [ ] Register insight with @register_insight
- [ ] Create optimization template JSON
- [ ] Write `SYS_17_DYNAMIC_RESPONSE.md` protocol
- [ ] Update `SYS_12_EXTRACTOR_LIBRARY.md`
- [ ] End-to-end test: optimization with random vib objective

### Phase 4+ (Future)
- [ ] Sine sweep processor
- [ ] Shock response spectrum
- [ ] Multi-axis SRSS/CQC
- [ ] Fatigue from random
- [ ] Neural surrogate for modal database

---

## 6. Success Criteria

### Technical
- [ ] Grms matches Miles' equation within 1% for SDOF
- [ ] Multi-mode response matches SOL 111 within 5%
- [ ] Processing time <5s per trial (excluding SOL 103)
- [ ] Memory <500MB for large models

### Integration
- [ ] `extract_random_response()` works in `run_optimization.py`
- [ ] DynamicResponseInsight appears in dashboard
- [ ] Protocol SYS_17 is complete and accurate
- [ ] Extractor library updated with E23-E25

### User Experience
- [ ] User can optimize for Grms using natural language
- [ ] Dashboard shows response PSD and Grms history
- [ ] Error messages are clear and actionable

---

## 7. Recommended Implementation Order

1. **Day 1-2**: Phase 1.1-1.4 (data classes, transfer functions, PSD)
2. **Day 3**: Phase 1.5 (tests) + Phase 2.1 (random processor core)
3. **Day 4**: Phase 2.2 (tests) + Phase 3.1 (extractor)
4. **Day 5**: Phase 3.2-3.4 (registry, insight, docs)
5. **Day 6**: End-to-end testing, validation, refinement

---

*This plan aligns the Master Plan with Atomizer's architecture for smooth implementation.*