Atomizer/.claude/skills/modules/DYNAMIC_RESPONSE_MASTER_PLAN.md

# Dynamic Response Processor - Master Implementation Plan

## Atomizer Feature: Modal Superposition Engine for Dynamic Optimization

**Document Version**: 1.0  
**Created**: 2025-12-22  
**Author**: Antoine + Claude  
**Status**: Implementation Ready  
**Priority**: High - Differentiating Feature

---

## Executive Summary

This document provides the complete implementation roadmap for adding **analytical dynamic response processing** to Atomizer. This feature enables optimization of structures under random vibration, sine sweep, and shock environments — **without running expensive frequency response FEA solutions for each design iteration**.

### Value Proposition

| Traditional Workflow | Atomizer Dynamic Workflow |
|---------------------|---------------------------|
| Run SOL 103 (modal) | Run SOL 103 (modal) |
| Run SOL 111 (frequency response) | **Skip** - compute analytically |
| Run SOL 112 (random) | **Skip** - compute analytically |
| ~15-30 min per iteration | ~30 sec per iteration |
| 100 iterations = 25-50 hours | 100 iterations = **50 minutes** |

**Key Differentiator**: SATK is a post-processor. Atomizer becomes an **optimization engine with built-in dynamic response intelligence**.

---

## Table of Contents

1. [Technical Foundation](#1-technical-foundation)
2. [Architecture Design](#2-architecture-design)
3. [Implementation Phases](#3-implementation-phases)
4. [Module Specifications](#4-module-specifications)
5. [Integration Points](#5-integration-points)
6. [Validation Strategy](#6-validation-strategy)
7. [Performance Targets](#7-performance-targets)
8. [Risk Mitigation](#8-risk-mitigation)
9. [Future Extensions](#9-future-extensions)

---

## 1. Technical Foundation

### 1.1 The Physics: Why This Works

The fundamental insight is that **dynamic response can be computed from modal data alone**.

For a structure with N degrees of freedom, the equation of motion is:

```
[M]{ü} + [C]{u̇} + [K]{u} = {F(t)}
```

Modal analysis (SOL 103) decomposes this into N independent single-DOF equations:

```
q̈ᵢ + 2ζᵢωᵢq̇ᵢ + ωᵢ²qᵢ = Lᵢ·a_base / mᵢ
```

Where:
- `qᵢ` = modal coordinate (how much mode i is excited)
- `ωᵢ` = natural frequency of mode i (rad/s)
- `ζᵢ` = damping ratio of mode i
- `Lᵢ` = modal participation factor (how mode i couples to base motion)
- `mᵢ` = generalized (modal) mass

**The key**: Once we extract (ωᵢ, ζᵢ, Lᵢ, φᵢ) from SOL 103, we can compute the response to ANY excitation profile analytically — no more FEA needed.

### 1.2 Transfer Function Mathematics

#### Single Mode Transfer Function

For base acceleration input, the displacement transfer function of mode i is:

```python
def H_displacement(f, f_n, zeta):
    """
    Displacement per unit base acceleration.
    
    H(f) = 1 / [(ωₙ² - ω²) + j·2·ζ·ωₙ·ω]
    
    In dimensionless form with r = f/fₙ:
    H(r) = 1 / [(1 - r²) + j·2·ζ·r]
    """
    r = f / f_n
    return 1 / ((1 - r**2) + 2j * zeta * r)
```

#### Acceleration Transmissibility

For base excitation, the ratio of response acceleration to input acceleration:

```python
def transmissibility(f, f_n, zeta):
    """
    T(f) = (1 + j·2·ζ·r) / [(1 - r²) + j·2·ζ·r]
    
    |T| = 1 at low frequency (rigid body)
    |T| = Q at resonance (amplification)
    |T| → 0 at high frequency (isolation)
    """
    r = f / f_n
    return (1 + 2j * zeta * r) / ((1 - r**2) + 2j * zeta * r)
```

#### Multi-Mode Combination

Total response sums all modal contributions:

```python
def total_response(frequencies, modes, participation_factors):
    """
    H_total(f) = Σᵢ Lᵢ · φᵢ · Hᵢ(f) / mᵢ
    
    Where:
    - Lᵢ = participation factor
    - φᵢ = mode shape at point of interest
    - Hᵢ = modal transfer function
    - mᵢ = generalized mass
    """
    H_total = np.zeros(len(frequencies), dtype=complex)
    for mode in modes:
        H_mode = H_displacement(frequencies, mode.frequency, mode.damping)
        H_total += mode.participation * mode.shape * H_mode / mode.gen_mass
    return H_total
```

### 1.3 Random Vibration Response

Random vibration is characterized by Power Spectral Density (PSD), typically in G²/Hz.

#### Response PSD

```python
S_response(f) = |H(f)|² × S_input(f)
```

#### RMS (Root Mean Square)

```python
σ² = ∫ S_response(f) df  # Mean square
RMS = √σ²                 # Root mean square
```

#### Peak Response

For Gaussian random processes:
```python
Peak = k × RMS
```

Where k (crest factor) depends on:
- **k = 3.0**: Standard 3-sigma (99.7% probability not exceeded)
- **k = 3.5-4.0**: Von Mises stress (non-Gaussian distribution)
- **k = 4.0-4.5**: Fatigue-critical applications

#### Miles' Equation (SDOF Approximation)

For quick estimates with white noise PSD:

```python
Grms = √(π/2 × fₙ × Q × PSD_level)
```

Where Q = 1/(2ζ) is the quality factor.

### 1.4 Stress Recovery from Modal Data

Element stresses can also be computed via modal superposition:

```python
σ(t) = Σᵢ qᵢ(t) × σᵢ_modal
```

Where `σᵢ_modal` is the stress pattern when mode i has unit amplitude.

For random vibration, stress PSD:
```python
S_σ(f) = |Σᵢ Lᵢ × σᵢ_modal × Hᵢ(f)|² × S_input(f)
```

### 1.5 What We Extract from SOL 103

| Data | Source | Used For |
|------|--------|----------|
| Natural frequencies (fₙ) | OP2: eigenvalues | Transfer function poles |
| Mode shapes (φ) | OP2: eigenvectors | Response at specific nodes |
| Generalized mass (m) | OP2/F06: eigenvalues | Modal scaling |
| Modal effective mass | F06: MEFFMASS | Participation estimation |
| Participation factors (L) | F06: MEFFMASS | Base excitation coupling |
| Modal stress shapes | OP2: modal stress | Stress recovery |

---

## 2. Architecture Design

### 2.1 Module Structure

```
optimization_engine/
├── processors/
│   └── dynamic_response/
│       ├── __init__.py                 # Public API exports
│       ├── modal_database.py           # ModalDatabase, ModeShape classes
│       ├── transfer_functions.py       # FRF computation engine
│       ├── random_vibration.py         # Random response processor
│       ├── sine_sweep.py               # Harmonic response processor
│       ├── shock_response.py           # SRS computation (Phase 5)
│       ├── stress_recovery.py          # Modal stress combination
│       ├── psd_profiles.py             # Standard PSD definitions
│       └── utils/
│           ├── __init__.py
│           ├── frequency_utils.py      # Frequency array generation
│           ├── integration.py          # Spectral integration methods
│           └── statistics.py           # Peak factors, distributions
│
├── extractors/
│   ├── extract_modal_database.py       # Full modal data extraction
│   ├── extract_random_response.py      # Random vibration objectives
│   ├── extract_sine_response.py        # Sine sweep objectives
│   └── extract_dynamic_stress.py       # Dynamic stress extraction
│
├── hooks/
│   └── dynamic_response/
│       ├── pre_dynamic_analysis.py     # Validate modal data
│       └── post_dynamic_analysis.py    # Cache results, logging
│
└── templates/
    ├── random_vibration_optimization.json
    ├── sine_sweep_optimization.json
    └── multi_environment_optimization.json
```

### 2.2 Class Hierarchy

```
┌─────────────────────────────────────────────────────────────┐
│                    DynamicResponseProcessor                  │
│                        (Abstract Base)                       │
├─────────────────────────────────────────────────────────────┤
│ + modal_db: ModalDatabase                                   │
│ + frf_engine: TransferFunctionEngine                        │
│ + compute_response() → ResponseResult                       │
│ + compute_at_node(node_id) → float                         │
│ + compute_stress(element_ids) → Dict                        │
└─────────────────────────────────────────────────────────────┘
                              │
          ┌───────────────────┼───────────────────┐
          ▼                   ▼                   ▼
┌─────────────────┐  ┌─────────────────┐  ┌─────────────────┐
│ RandomVibration │  │   SineSweep     │  │ ShockResponse   │
│   Processor     │  │   Processor     │  │   Processor     │
├─────────────────┤  ├─────────────────┤  ├─────────────────┤
│ + psd_profile   │  │ + sweep_rate    │  │ + srs_damping   │
│ + compute_grms()│  │ + dwell_freqs   │  │ + compute_srs() │
│ + compute_peak()│  │ + compute_peak()│  │ + compute_pseudo│
└─────────────────┘  └─────────────────┘  └─────────────────┘
```

### 2.3 Data Flow

```
┌──────────────────────────────────────────────────────────────────┐
│                      OPTIMIZATION LOOP                            │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌──────────────────────────────────────────────────────────────────┐
│  1. UPDATE DESIGN PARAMETERS                                      │
│     NX Expression Manager → thickness=2.5mm, rib_height=12mm     │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌──────────────────────────────────────────────────────────────────┐
│  2. RUN SOL 103 (Modal Analysis)                                  │
│     NX Solver → model.op2, model.f06                             │
│     Time: ~20-60 seconds                                          │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌──────────────────────────────────────────────────────────────────┐
│  3. EXTRACT MODAL DATABASE                                        │
│     pyNastran → ModalDatabase object                              │
│     - Frequencies, mode shapes, participation factors             │
│     - Modal stress shapes (if needed)                             │
│     Time: ~0.5-2 seconds                                          │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌──────────────────────────────────────────────────────────────────┐
│  4. COMPUTE DYNAMIC RESPONSE (Analytical)                         │
│     RandomVibrationProcessor.compute_response()                   │
│     - Build transfer functions: ~5ms                              │
│     - Compute response PSD: ~10ms                                 │
│     - Integrate for RMS: ~5ms                                     │
│     - Stress recovery: ~50ms                                      │
│     Time: ~50-100 milliseconds                                    │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌──────────────────────────────────────────────────────────────────┐
│  5. RETURN OBJECTIVES & CONSTRAINTS                               │
│     {                                                             │
│       'grms': 8.4,                                               │
│       'peak_acceleration': 25.2,                                  │
│       'peak_stress': 145.6,                                       │
│       'first_frequency': 127.3                                    │
│     }                                                             │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌──────────────────────────────────────────────────────────────────┐
│  6. OPTUNA SUGGESTS NEXT DESIGN                                   │
│     TPE/CMA-ES/NSGA-II → next trial parameters                   │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
                        [REPEAT 1-6]
```

### 2.4 Memory Management

For large models with many modes/elements, we need smart memory handling:

```python
class ModalDatabase:
    """
    Memory-efficient modal data storage.
    
    Strategies:
    1. Lazy loading - only load mode shapes when needed
    2. Sparse storage - only store non-zero eigenvector components
    3. Node filtering - only extract nodes of interest
    4. Disk caching - cache to HDF5 for reuse across trials
    """
    
    def __init__(self, op2_file, lazy_load=True):
        self._op2_path = op2_file
        self._modes_loaded = False
        self._mode_shapes = None  # Load on demand
        
    @property
    def mode_shapes(self):
        if not self._modes_loaded:
            self._load_mode_shapes()
        return self._mode_shapes
    
    def get_modes_for_nodes(self, node_ids: List[int]) -> np.ndarray:
        """Extract mode shapes only for specified nodes."""
        # Efficient partial loading
        pass
    
    def cache_to_hdf5(self, cache_path: Path):
        """Cache modal database for fast reload."""
        pass
    
    @classmethod
    def from_cache(cls, cache_path: Path) -> 'ModalDatabase':
        """Load from HDF5 cache (10x faster than OP2)."""
        pass
```

---

## 3. Implementation Phases

### Phase 1: Core Infrastructure (Week 1)

**Goal**: Build the foundational classes and basic functionality.

#### Tasks

| Task | Description | Effort | Deliverable |
|------|-------------|--------|-------------|
| 1.1 | Create `ModeShape` dataclass | 2h | `modal_database.py` |
| 1.2 | Create `ModalDatabase` class | 4h | `modal_database.py` |
| 1.3 | Implement OP2 extraction | 4h | `ModalDatabase.from_op2_f06()` |
| 1.4 | Implement F06 MEFFMASS parsing | 3h | Modal effective mass extraction |
| 1.5 | Create `TransferFunctionEngine` | 4h | `transfer_functions.py` |
| 1.6 | Single mode FRF | 2h | `single_mode_frf()` |
| 1.7 | Multi-mode combination | 2h | `base_excitation_frf()` |
| 1.8 | Unit tests for FRF | 3h | `tests/test_transfer_functions.py` |

**Validation Checkpoint**: 
- FRF matches closed-form SDOF solution
- Modal database correctly extracts from sample OP2

#### Code: ModeShape Dataclass

```python
# optimization_engine/processors/dynamic_response/modal_database.py

from dataclasses import dataclass, field
from typing import Dict, Optional, List
import numpy as np

@dataclass
class ModeShape:
    """
    Single mode shape data container.
    
    Stores all information needed for modal superposition
    of one natural mode.
    """
    # Identification
    mode_number: int
    
    # Dynamic properties
    frequency: float  # Hz
    damping_ratio: float = 0.02  # Dimensionless (2% default)
    
    # Modal mass/stiffness
    generalized_mass: float = 1.0  # kg (normalized)
    generalized_stiffness: float = 0.0  # N/m
    
    # Participation factors (for base excitation)
    participation_x: float = 0.0
    participation_y: float = 0.0
    participation_z: float = 0.0
    
    # Modal effective mass (absolute contribution)
    effective_mass_x: float = 0.0
    effective_mass_y: float = 0.0
    effective_mass_z: float = 0.0
    
    # Eigenvector at specific nodes
    # Dict[node_id, np.array([tx, ty, tz, rx, ry, rz])]
    eigenvector: Dict[int, np.ndarray] = field(default_factory=dict)
    
    # Modal stress shapes 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)
    
    # Derived properties
    @property
    def omega(self) -> float:
        """Angular frequency (rad/s)."""
        return 2 * np.pi * self.frequency
    
    @property
    def omega_squared(self) -> float:
        """ω² for efficiency."""
        return self.omega ** 2
    
    @property
    def quality_factor(self) -> float:
        """Q = 1/(2ζ) - amplification at resonance."""
        if self.damping_ratio > 0:
            return 1 / (2 * self.damping_ratio)
        return np.inf
    
    @property
    def half_power_bandwidth(self) -> float:
        """Δf = fₙ/Q - bandwidth at -3dB."""
        return self.frequency / self.quality_factor
    
    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:
        """Get effective mass for direction."""
        return getattr(self, f'effective_mass_{direction}', 0.0)
    
    def get_eigenvector_component(
        self, 
        node_id: int, 
        direction: str
    ) -> float:
        """Get eigenvector component at node in direction."""
        if node_id not in self.eigenvector:
            return 0.0
        dir_map = {'x': 0, 'y': 1, 'z': 2, 'rx': 3, 'ry': 4, 'rz': 5}
        idx = dir_map.get(direction, 0)
        return self.eigenvector[node_id][idx]
```

### Phase 2: Random Vibration Processor (Week 2)

**Goal**: Full random vibration response calculation with stress recovery.

#### Tasks

| Task | Description | Effort | Deliverable |
|------|-------------|--------|-------------|
| 2.1 | Create `PSDProfile` class | 2h | `psd_profiles.py` |
| 2.2 | Standard profiles (GEVS, MIL-STD) | 2h | `STANDARD_PROFILES` dict |
| 2.3 | PSD interpolation (log-log) | 2h | `PSDProfile.interpolate()` |
| 2.4 | Create `RandomVibrationProcessor` | 4h | `random_vibration.py` |
| 2.5 | Response PSD calculation | 3h | `compute_response_psd()` |
| 2.6 | Spectral integration | 2h | `compute_rms()` |
| 2.7 | Peak factor calculation | 2h | Gaussian + non-Gaussian |
| 2.8 | Modal contribution analysis | 2h | Identify dominant modes |
| 2.9 | Stress recovery | 4h | `compute_stress_response()` |
| 2.10 | Result dataclass | 2h | `RandomVibrationResult` |
| 2.11 | Integration tests | 4h | Miles' equation validation |

**Validation Checkpoint**:
- Grms matches Miles' equation for SDOF
- Multi-mode response reasonable for sample model
- Stress recovery gives correct distribution

#### Code: PSD Profiles

```python
# optimization_engine/processors/dynamic_response/psd_profiles.py

import numpy as np
from scipy import interpolate
from dataclasses import dataclass
from typing import List, Tuple, Optional
from pathlib import Path


@dataclass
class PSDProfile:
    """
    Power Spectral Density profile for random vibration.
    
    Stores PSD as G²/Hz vs frequency, with log-log interpolation.
    """
    name: str
    frequencies: np.ndarray  # Hz
    psd_values: np.ndarray   # G²/Hz
    description: str = ""
    source: str = ""  # Standard reference
    
    @property
    def freq_min(self) -> float:
        return float(self.frequencies.min())
    
    @property
    def freq_max(self) -> float:
        return float(self.frequencies.max())
    
    @property
    def grms(self) -> float:
        """Overall Grms level."""
        from scipy.integrate import trapezoid
        return np.sqrt(trapezoid(self.psd_values, self.frequencies))
    
    def interpolate(self, frequencies: np.ndarray) -> np.ndarray:
        """
        Interpolate PSD to arbitrary frequencies.
        
        Uses log-log interpolation (standard for PSD).
        Clamps to edge values outside defined range.
        """
        # Safety: ensure positive frequencies
        frequencies = np.maximum(frequencies, 1e-6)
        
        log_f = np.log10(self.frequencies)
        log_psd = np.log10(np.maximum(self.psd_values, 1e-20))
        
        interp_func = interpolate.interp1d(
            log_f, log_psd,
            kind='linear',
            bounds_error=False,
            fill_value=(log_psd[0], log_psd[-1])
        )
        
        return 10 ** interp_func(np.log10(frequencies))
    
    def plot(self, ax=None, **kwargs):
        """Plot PSD profile."""
        import matplotlib.pyplot as plt
        if ax is None:
            fig, ax = plt.subplots()
        ax.loglog(self.frequencies, self.psd_values, **kwargs)
        ax.set_xlabel('Frequency (Hz)')
        ax.set_ylabel('PSD (G²/Hz)')
        ax.set_title(f'{self.name} - {self.grms:.2f} Grms')
        ax.grid(True, which='both', alpha=0.3)
        return ax
    
    @classmethod
    def from_breakpoints(
        cls, 
        breakpoints: List[Tuple[float, float]],
        name: str = "Custom",
        description: str = ""
    ) -> 'PSDProfile':
        """
        Create PSD from breakpoint specification.
        
        Args:
            breakpoints: List of (frequency_hz, psd_g2_hz) tuples
            name: Profile name
            description: Description text
        
        Example:
            psd = PSDProfile.from_breakpoints([
                (20, 0.01),    # 20 Hz: 0.01 G²/Hz
                (80, 0.04),    # 80 Hz: 0.04 G²/Hz (ramp up)
                (500, 0.04),   # 500 Hz: 0.04 G²/Hz (flat)
                (2000, 0.007)  # 2000 Hz: 0.007 G²/Hz (ramp down)
            ], name='MIL-STD-810')
        """
        breakpoints = sorted(breakpoints, key=lambda x: x[0])
        frequencies = np.array([bp[0] for bp in breakpoints])
        psd_values = np.array([bp[1] for bp in breakpoints])
        return cls(
            name=name,
            frequencies=frequencies,
            psd_values=psd_values,
            description=description
        )
    
    @classmethod
    def from_csv(cls, filepath: Path, name: str = None) -> 'PSDProfile':
        """Load PSD from CSV file with freq,psd columns."""
        data = np.loadtxt(filepath, delimiter=',', skiprows=1)
        return cls(
            name=name or filepath.stem,
            frequencies=data[:, 0],
            psd_values=data[:, 1]
        )
    
    @classmethod
    def flat(
        cls,
        level: float,
        freq_min: float = 20,
        freq_max: float = 2000,
        name: str = "Flat"
    ) -> 'PSDProfile':
        """Create flat (white noise) PSD profile."""
        return cls(
            name=name,
            frequencies=np.array([freq_min, freq_max]),
            psd_values=np.array([level, level]),
            description=f"Flat PSD at {level} G²/Hz"
        )


# ============================================================================
# STANDARD PSD PROFILES
# ============================================================================

STANDARD_PROFILES = {
    # NASA GEVS (General Environmental Verification Standard)
    'GEVS': PSDProfile.from_breakpoints([
        (20, 0.026),
        (50, 0.16),
        (800, 0.16),
        (2000, 0.026),
    ], name='NASA GEVS', 
       description='NASA General Environmental Verification Standard'),
    
    # MIL-STD-810G Category 1 (Basic)
    'MIL-STD-810G_CAT1': PSDProfile.from_breakpoints([
        (10, 0.04),
        (40, 0.04),
        (500, 0.04),
        (2000, 0.01),
    ], name='MIL-STD-810G Cat 1',
       description='Basic transportation vibration'),
    
    # MIL-STD-810G Category 4 (Truck)
    'MIL-STD-810G_CAT4': PSDProfile.from_breakpoints([
        (5, 0.001),
        (10, 0.015),
        (40, 0.015),
        (500, 0.015),
        (2000, 0.003),
    ], name='MIL-STD-810G Cat 4',
       description='Truck transportation'),
    
    # MIL-STD-810G Category 7 (Jet Aircraft)
    'MIL-STD-810G_CAT7': PSDProfile.from_breakpoints([
        (15, 0.015),
        (80, 0.015),
        (350, 0.015),
        (2000, 0.006),
    ], name='MIL-STD-810G Cat 7',
       description='Jet aircraft vibration'),
    
    # MIL-STD-810G Category 24 (Helicopter)
    'MIL-STD-810G_CAT24': PSDProfile.from_breakpoints([
        (5, 0.001),
        (20, 0.01),
        (80, 0.04),
        (350, 0.04),
        (2000, 0.007),
    ], name='MIL-STD-810G Cat 24',
       description='Helicopter vibration'),
    
    # NAVMAT P-9492
    'NAVMAT': PSDProfile.from_breakpoints([
        (20, 0.04),
        (80, 0.04),
        (350, 0.04),
        (2000, 0.007),
    ], name='NAVMAT P-9492',
       description='Navy shipboard equipment'),
    
    # ESA PSS-01-401 Qualification
    'ESA_QUAL': PSDProfile.from_breakpoints([
        (20, 0.01),
        (100, 0.08),
        (300, 0.08),
        (2000, 0.008),
    ], name='ESA PSS-01-401',
       description='ESA qualification level'),
    
    # ESA PSS-01-401 Acceptance
    'ESA_ACCEPT': PSDProfile.from_breakpoints([
        (20, 0.005),
        (100, 0.04),
        (300, 0.04),
        (2000, 0.004),
    ], name='ESA PSS-01-401 Accept',
       description='ESA acceptance level'),
    
    # SpaceX Falcon 9 (approximate public data)
    'FALCON9': PSDProfile.from_breakpoints([
        (20, 0.01),
        (50, 0.08),
        (800, 0.08),
        (2000, 0.02),
    ], name='Falcon 9 (approx)',
       description='SpaceX Falcon 9 approximate envelope'),
    
    # Ariane 5 (approximate)
    'ARIANE5': PSDProfile.from_breakpoints([
        (20, 0.015),
        (100, 0.09),
        (400, 0.09),
        (2000, 0.01),
    ], name='Ariane 5 (approx)',
       description='Ariane 5 approximate envelope'),
}


def get_psd_profile(name: str) -> PSDProfile:
    """
    Get a standard PSD profile by name.
    
    Args:
        name: Profile name (case-insensitive, underscores/hyphens flexible)
    
    Returns:
        PSDProfile instance
    
    Raises:
        ValueError: If profile not found
    
    Example:
        >>> psd = get_psd_profile('GEVS')
        >>> psd = get_psd_profile('mil-std-810g-cat7')
    """
    # Normalize name
    normalized = name.upper().replace('-', '_').replace(' ', '_')
    
    if normalized in STANDARD_PROFILES:
        return STANDARD_PROFILES[normalized]
    
    # Try partial match
    for key, profile in STANDARD_PROFILES.items():
        if normalized in key or key in normalized:
            return profile
    
    available = ', '.join(STANDARD_PROFILES.keys())
    raise ValueError(f"Unknown PSD profile: '{name}'. Available: {available}")


def list_available_profiles() -> List[str]:
    """Return list of available standard PSD profile names."""
    return list(STANDARD_PROFILES.keys())
```

### Phase 3: Atomizer Integration (Week 3)

**Goal**: Connect dynamic response processor to Atomizer's extractor and optimization systems.

#### Tasks

| Task | Description | Effort | Deliverable |
|------|-------------|--------|-------------|
| 3.1 | Create `extract_modal_database` | 3h | New extractor |
| 3.2 | Create `extract_random_response` | 3h | New extractor |
| 3.3 | Create `extract_dynamic_stress` | 3h | New extractor |
| 3.4 | Optimization template | 2h | `random_vibration_optimization.json` |
| 3.5 | Runner integration | 4h | Hook into `OptimizationRunner` |
| 3.6 | Dashboard widgets | 4h | Response PSD plot, Grms history |
| 3.7 | Result caching | 3h | Cache modal DB between similar designs |
| 3.8 | Documentation | 3h | User guide, API docs |
| 3.9 | End-to-end test | 4h | Full optimization run |

**Validation Checkpoint**:
- Complete optimization runs successfully
- Dashboard shows dynamic response metrics
- Results match manual SATK/Nastran verification

#### Code: Random Response Extractor

```python
# optimization_engine/extractors/extract_random_response.py

"""
Random Vibration Response Extractor for Atomizer Optimization
=============================================================

Extracts random vibration metrics (Grms, peak acceleration, peak stress)
from SOL 103 modal analysis results using analytical modal superposition.

This extractor enables dynamic optimization without running SOL 111.

Example:
    >>> result = extract_random_response(
    ...     op2_file='model.op2',
    ...     psd_profile='GEVS',
    ...     direction='z'
    ... )
    >>> print(f"Grms: {result['grms']:.2f} G")
    >>> print(f"Peak: {result['peak_acceleration']:.2f} G")
"""

from pathlib import Path
from typing import Dict, Any, Optional, Union, List
import logging

logger = logging.getLogger(__name__)


def extract_random_response(
    op2_file: Union[str, Path],
    psd_profile: str = 'GEVS',
    direction: str = 'z',
    damping: float = 0.02,
    node_id: Optional[int] = None,
    peak_sigma: float = 3.0,
    f06_file: Optional[Union[str, Path]] = None,
    n_modes: Optional[int] = None,
    frequency_resolution: int = 500,
) -> Dict[str, Any]:
    """
    Extract random vibration response for optimization objective/constraint.
    
    This is the main interface for dynamic optimization. It reads SOL 103
    results and computes random vibration response analytically.
    
    Args:
        op2_file: Path to OP2 file from SOL 103 analysis
        psd_profile: PSD specification - name ('GEVS', 'MIL-STD-810G_CAT7', etc.)
                     or path to CSV file
        direction: Excitation direction ('x', 'y', 'z')
        damping: Modal damping ratio (default 2%)
        node_id: Specific node for response (None = use effective mass)
        peak_sigma: Sigma level for peak (3.0 = 3-sigma)
        f06_file: Optional F06 for modal effective mass table
        n_modes: Number of modes to include (None = auto 95% mass)
        frequency_resolution: Points in frequency array
    
    Returns:
        Dict containing:
        - 'grms': RMS acceleration (G)
        - 'peak_acceleration': Peak acceleration (G)
        - 'rms_displacement': RMS displacement (mm)
        - 'peak_displacement': Peak displacement (mm)
        - 'dominant_mode': Mode number with highest contribution
        - 'dominant_frequency': Frequency of dominant mode (Hz)
        - 'n_modes_used': Number of modes in calculation
        - 'mass_fraction': Cumulative mass fraction of included modes
        - 'psd_profile': Name of PSD profile used
        - 'direction': Excitation direction
        - 'success': True if extraction succeeded
        - 'error': Error message if failed
    
    Example:
        # For optimization objective
        >>> result = extract_random_response('bracket.op2', psd_profile='GEVS')
        >>> objective = result['peak_acceleration']
        
        # For constraint checking
        >>> result = extract_random_response('bracket.op2', psd_profile='GEVS')
        >>> if result['peak_acceleration'] > 25.0:
        ...     print("Constraint violated!")
    """
    from optimization_engine.processors.dynamic_response import (
        ModalDatabase,
        RandomVibrationProcessor,
        get_psd_profile
    )
    
    op2_path = Path(op2_file)
    f06_path = Path(f06_file) if f06_file else None
    
    # Auto-detect F06 if not provided
    if f06_path is None:
        possible_f06 = op2_path.with_suffix('.f06')
        if possible_f06.exists():
            f06_path = possible_f06
            logger.info(f"Auto-detected F06: {f06_path}")
    
    try:
        # Build modal database
        logger.info(f"Extracting modal database from {op2_path}")
        modal_db = ModalDatabase.from_op2_f06(
            op2_file=op2_path,
            f06_file=f06_path,
            damping=damping,
            node_ids=[node_id] if node_id else None,
        )
        
        logger.info(f"Extracted {modal_db.n_modes} modes, "
                   f"freq range: {modal_db.frequency_range[0]:.1f}-"
                   f"{modal_db.frequency_range[1]:.1f} Hz")
        
        # Load PSD profile
        if Path(psd_profile).exists():
            from optimization_engine.processors.dynamic_response import PSDProfile
            psd = PSDProfile.from_csv(Path(psd_profile))
        else:
            psd = get_psd_profile(psd_profile)
        
        logger.info(f"Using PSD profile: {psd.name} ({psd.grms:.2f} Grms)")
        
        # Create processor
        processor = RandomVibrationProcessor(
            modal_db=modal_db,
            default_damping=damping
        )
        
        # Compute response
        result = processor.compute_response(
            psd_profile=psd,
            direction=direction,
            node_id=node_id,
            n_modes=n_modes,
            frequency_resolution=frequency_resolution,
            peak_sigma=peak_sigma,
            output_type='acceleration'
        )
        
        # Get mass fraction for included modes
        n_modes_used = n_modes or modal_db.modes_for_mass_fraction(0.95, direction)
        mass_fraction = modal_db.get_cumulative_mass_fraction(direction)[n_modes_used - 1]
        
        # Dominant mode info
        dominant_freq = modal_db.modes[result.dominant_mode - 1].frequency
        
        return {
            'grms': result.rms_acceleration,
            'peak_acceleration': result.peak_acceleration,
            'rms_displacement': result.rms_displacement,
            'peak_displacement': result.peak_displacement,
            'dominant_mode': result.dominant_mode,
            'dominant_frequency': dominant_freq,
            'n_modes_used': n_modes_used,
            'mass_fraction': mass_fraction,
            'modal_contributions': result.modal_contributions,
            'psd_profile': psd.name,
            'psd_grms': psd.grms,
            'direction': direction,
            'peak_sigma': peak_sigma,
            'success': True,
            'error': None,
        }
        
    except Exception as e:
        logger.exception(f"Random response extraction failed: {e}")
        return {
            'grms': None,
            'peak_acceleration': None,
            'rms_displacement': None,
            'peak_displacement': None,
            'dominant_mode': None,
            'dominant_frequency': None,
            'success': False,
            'error': str(e),
        }


def extract_random_stress(
    op2_file: Union[str, Path],
    psd_profile: str = 'GEVS',
    direction: str = 'z',
    damping: float = 0.02,
    peak_sigma: float = 3.0,
    element_ids: Optional[List[int]] = None,
    non_gaussian_correction: bool = True,
    f06_file: Optional[Union[str, Path]] = None,
) -> Dict[str, Any]:
    """
    Extract peak Von Mises stress from random vibration.
    
    Args:
        op2_file: Path to OP2 with modal stress results
        psd_profile: PSD specification
        direction: Excitation direction
        damping: Modal damping ratio
        peak_sigma: Sigma level
        element_ids: Specific elements (None = all)
        non_gaussian_correction: Apply VM stress distribution correction
        f06_file: Optional F06 for effective mass
    
    Returns:
        Dict with:
        - 'max_rms_von_mises': Maximum RMS VM stress (MPa)
        - 'max_peak_von_mises': Maximum peak VM stress (MPa)
        - 'critical_element': Element ID with max stress
        - 'margin_of_safety': (allowable - peak) / peak
        - 'n_elements': Number of elements analyzed
    """
    from optimization_engine.processors.dynamic_response import (
        ModalDatabase,
        RandomVibrationProcessor,
        get_psd_profile
    )
    
    op2_path = Path(op2_file)
    f06_path = Path(f06_file) if f06_file else None
    
    try:
        # Build modal database with stress shapes
        modal_db = ModalDatabase.from_op2_f06(
            op2_file=op2_path,
            f06_file=f06_path,
            damping=damping,
            element_ids=element_ids,  # Request stress shapes
        )
        
        # Load PSD
        psd = get_psd_profile(psd_profile) if not Path(psd_profile).exists() \
              else PSDProfile.from_csv(Path(psd_profile))
        
        # Process
        processor = RandomVibrationProcessor(modal_db, damping)
        result = processor.compute_stress_response(
            psd_profile=psd,
            direction=direction,
            element_ids=element_ids,
            peak_sigma=peak_sigma,
            non_gaussian_correction=non_gaussian_correction,
        )
        
        return {
            'max_rms_von_mises': result['max_rms_von_mises'],
            'max_peak_von_mises': result['max_peak_von_mises'],
            'critical_element': result['critical_element'],
            'n_elements': result['n_elements'],
            'success': True,
            'error': None,
        }
        
    except Exception as e:
        logger.exception(f"Stress extraction failed: {e}")
        return {
            'max_rms_von_mises': None,
            'max_peak_von_mises': None,
            'critical_element': None,
            'success': False,
            'error': str(e),
        }


# Convenience functions for optimization objectives

def get_grms(op2_file, psd='GEVS', direction='z', damping=0.02) -> float:
    """Quick Grms extraction for objective function."""
    result = extract_random_response(op2_file, psd, direction, damping)
    return result['grms'] if result['success'] else float('inf')


def get_peak_g(op2_file, psd='GEVS', direction='z', damping=0.02, sigma=3.0) -> float:
    """Quick peak acceleration extraction for constraint."""
    result = extract_random_response(op2_file, psd, direction, damping, peak_sigma=sigma)
    return result['peak_acceleration'] if result['success'] else float('inf')


def get_peak_stress(op2_file, psd='GEVS', direction='z', damping=0.02, sigma=3.0) -> float:
    """Quick peak stress extraction for constraint."""
    result = extract_random_stress(op2_file, psd, direction, damping, sigma)
    return result['max_peak_von_mises'] if result['success'] else float('inf')
```

### Phase 4: Sine Sweep Processor (Week 4)

**Goal**: Add harmonic/sine sweep response capability.

#### Tasks

| Task | Description | Effort | Deliverable |
|------|-------------|--------|-------------|
| 4.1 | Sine sweep transfer function | 3h | Magnitude/phase at frequencies |
| 4.2 | Resonance search | 2h | Find peak frequencies |
| 4.3 | Dwell analysis | 3h | Response at specific frequencies |
| 4.4 | Sweep rate effects | 2h | Transient buildup correction |
| 4.5 | Create `SineSweepProcessor` | 4h | `sine_sweep.py` |
| 4.6 | Create extractor | 2h | `extract_sine_response.py` |
| 4.7 | Optimization template | 2h | `sine_sweep_optimization.json` |
| 4.8 | Tests | 3h | Validate against analytical |

### Phase 5: Advanced Features (Future)

**Goal**: Extended capabilities for complete dynamic qualification.

#### Shock Response Spectrum (SRS)
- Compute SRS from transient pulse or base motion
- Pseudo-velocity and pseudo-acceleration
- Compare against equipment SRS limits

#### Multi-Axis Excitation
- SRSS (Square Root Sum of Squares) combination
- CQC (Complete Quadratic Combination) for correlated inputs
- Full 6-DOF base motion

#### Fatigue from Random Vibration
- Narrow-band fatigue using Miner's rule
- Dirlik method for wide-band random
- Cycle counting from PSD

#### Neural Surrogate for Modal Database
- Train GNN to predict (frequencies, mode shapes) from geometry
- Skip SOL 103 entirely for similar designs
- 1000x additional speedup potential

---

## 4. Module Specifications

### 4.1 ModalDatabase

```python
class ModalDatabase:
    """
    Central storage for all modal analysis data.
    
    Responsibilities:
    - Extract data from OP2/F06 files
    - Store mode shapes, frequencies, participation factors
    - Provide efficient access patterns for dynamic calculations
    - Support caching and serialization
    
    Thread Safety: Read-only after construction
    Memory: O(n_modes × n_nodes × 6) for full eigenvectors
    """
    
    # Class attributes
    SUPPORTED_ELEMENTS = ['CQUAD4', 'CTRIA3', 'CHEXA', 'CPENTA', 'CTETRA']
    
    # Required data
    model_name: str
    total_mass: float
    center_of_gravity: np.ndarray
    modes: List[ModeShape]
    
    # Optional data
    node_coordinates: Dict[int, np.ndarray]
    element_connectivity: Dict[int, List[int]]
    
    # Metadata
    source_op2: Path
    source_f06: Path
    extraction_timestamp: datetime
    nastran_version: str
    
    # Methods
    @classmethod
    def from_op2_f06(cls, op2_file, f06_file, **options) -> 'ModalDatabase'
    
    def get_mode(self, mode_number: int) -> ModeShape
    def get_frequencies(self) -> np.ndarray
    def get_participation_factors(self, direction: str) -> np.ndarray
    def get_cumulative_mass_fraction(self, direction: str) -> np.ndarray
    def modes_for_mass_fraction(self, fraction: float, direction: str) -> int
    
    def to_hdf5(self, filepath: Path)
    @classmethod
    def from_hdf5(cls, filepath: Path) -> 'ModalDatabase'
    
    def to_dict(self) -> Dict
    def summary(self) -> str
```

### 4.2 TransferFunctionEngine

```python
class TransferFunctionEngine:
    """
    Computes frequency response functions from modal data.
    
    Responsibilities:
    - Single mode FRF calculation
    - Multi-mode combination
    - Base excitation transmissibility
    - Point excitation mobility (future)
    
    Performance: Vectorized numpy operations
    Accuracy: Double precision complex arithmetic
    """
    
    def __init__(self, modal_db: ModalDatabase)
    
    # Core FRF methods
    def single_mode_frf(
        self,
        frequencies: np.ndarray,
        mode: ModeShape,
        output_type: str = 'displacement'
    ) -> np.ndarray  # Complex FRF
    
    def base_excitation_frf(
        self,
        frequencies: np.ndarray,
        direction: str,
        output_type: str = 'acceleration',
        node_id: int = None,
        n_modes: int = None
    ) -> np.ndarray  # Complex transmissibility
    
    def point_excitation_frf(
        self,
        frequencies: np.ndarray,
        input_node: int,
        input_direction: str,
        output_node: int,
        output_direction: str
    ) -> np.ndarray  # Complex FRF
    
    # Stress FRF
    def stress_frf(
        self,
        frequencies: np.ndarray,
        element_id: int,
        direction: str
    ) -> np.ndarray  # Complex stress transfer function
```

### 4.3 RandomVibrationProcessor

```python
class RandomVibrationProcessor:
    """
    Computes random vibration response using modal superposition.
    
    This is the SATK-equivalent processor.
    
    Capabilities:
    - Grms and peak response calculation
    - Response PSD computation
    - Modal contribution analysis
    - Stress recovery with non-Gaussian correction
    - Multiple response locations
    
    Performance Target: <100ms for typical aerospace model
    """
    
    def __init__(self, modal_db: ModalDatabase, default_damping: float = 0.02)
    
    # Main computation methods
    def compute_response(
        self,
        psd_profile: PSDProfile,
        direction: str = 'z',
        node_id: int = None,
        n_modes: int = None,
        frequency_resolution: int = 500,
        peak_sigma: float = 3.0,
        output_type: str = 'acceleration'
    ) -> RandomVibrationResult
    
    def compute_stress_response(
        self,
        psd_profile: PSDProfile,
        direction: str = 'z',
        element_ids: List[int] = None,
        peak_sigma: float = 3.0,
        non_gaussian_correction: bool = True
    ) -> Dict[str, Any]
    
    # Quick objective/constraint methods
    def compute_grms(self, psd_profile, direction, node_id=None) -> float
    def compute_peak_g(self, psd_profile, direction, sigma=3.0) -> float
    def compute_peak_stress(self, psd_profile, direction, sigma=3.0) -> float
    
    # Analysis methods
    def get_modal_contributions(self, psd_profile, direction) -> Dict[int, float]
    def get_response_psd(self, psd_profile, direction) -> Tuple[np.ndarray, np.ndarray]
    def identify_critical_modes(self, psd_profile, direction, threshold=0.1) -> List[int]
```

---

## 5. Integration Points

### 5.1 Extractor Registry

Add to `optimization_engine/extractors/__init__.py`:

```python
# Dynamic response extractors
from .extract_modal_database import extract_modal_database
from .extract_random_response import (
    extract_random_response,
    extract_random_stress,
    get_grms,
    get_peak_g,
    get_peak_stress
)
from .extract_sine_response import extract_sine_response

__all__ = [
    # ... existing extractors ...
    
    # Dynamic response
    'extract_modal_database',
    'extract_random_response',
    'extract_random_stress',
    'get_grms',
    'get_peak_g', 
    'get_peak_stress',
    'extract_sine_response',
]
```

### 5.2 Feature Registry

Add to `optimization_engine/feature_registry.json`:

```json
{
  "dynamic_response": {
    "random_vibration_processor": {
      "feature_id": "random_vibration_processor",
      "name": "Random Vibration Processor",
      "description": "Modal superposition-based random vibration response calculation",
      "category": "dynamics",
      "implementation": {
        "file_path": "optimization_engine/processors/dynamic_response/random_vibration.py",
        "class_name": "RandomVibrationProcessor"
      },
      "interface": {
        "inputs": [
          {"name": "modal_database", "type": "ModalDatabase"},
          {"name": "psd_profile", "type": "PSDProfile"},
          {"name": "direction", "type": "str"}
        ],
        "outputs": [
          {"name": "grms", "type": "float"},
          {"name": "peak_acceleration", "type": "float"},
          {"name": "response_psd", "type": "np.ndarray"}
        ]
      }
    }
  }
}
```

### 5.3 Dashboard Integration

New dashboard components for dynamic response:

```typescript
// atomizer-dashboard/src/components/DynamicResponse/
├── ResponsePSDPlot.tsx      // Log-log PSD plot with input/output
├── TransmissibilityPlot.tsx // FRF magnitude vs frequency
├── ModalContributionBar.tsx // Bar chart of mode contributions
├── GrmsHistoryChart.tsx     // Grms vs trial number
├── CampbellDiagram.tsx      // For rotating machinery (future)
└── DynamicResponseCard.tsx  // Summary card for dashboard
```

### 5.4 CLI Integration

```bash
# New CLI commands
atomizer dynamic extract model.op2 --psd GEVS --direction z
atomizer dynamic psd-list
atomizer dynamic psd-plot GEVS --output psd.png
atomizer optimize --template random_vibration_optimization.json
```

---

## 6. Validation Strategy

### 6.1 Unit Tests

```python
# tests/test_dynamic_response/

test_modal_database.py
├── test_mode_shape_creation
├── test_op2_extraction
├── test_f06_parsing
├── test_mass_fraction_calculation
└── test_hdf5_caching

test_transfer_functions.py
├── test_sdof_frf_analytical
├── test_resonance_amplification
├── test_phase_at_resonance
├── test_multi_mode_combination
└── test_base_excitation_transmissibility

test_random_vibration.py
├── test_miles_equation_sdof
├── test_grms_calculation
├── test_peak_factor
├── test_non_gaussian_correction
├── test_modal_contributions
└── test_stress_recovery

test_psd_profiles.py
├── test_interpolation
├── test_grms_integration
├── test_standard_profiles
└── test_csv_loading
```

### 6.2 Integration Tests

```python
# tests/integration/

test_dynamic_optimization.py
├── test_full_random_vib_optimization
├── test_multi_objective_dynamic
├── test_constraint_handling
└── test_dashboard_updates

test_nastran_comparison.py
├── test_vs_sol111_response
├── test_vs_sol112_random
└── test_stress_recovery_accuracy
```

### 6.3 Validation Cases

| Case | Model | Validation Against | Acceptance |
|------|-------|-------------------|------------|
| SDOF | Single mass-spring | Miles' equation | <1% error |
| Cantilever beam | Simple beam | SOL 111 FRF | <5% error |
| Plate | Rectangular plate | SOL 112 random | <5% error |
| Bracket | Real component | SATK results | <10% error |

### 6.4 Benchmark Against SATK

If possible, run same model through both Atomizer and SATK:

```
Model: CubeSat bracket
PSD: GEVS
Direction: Z-axis
Damping: 2%

Compare:
- Grms
- Peak acceleration
- Dominant mode
- Response PSD shape
- Peak stress location
```

---

## 7. Performance Targets

### 7.1 Timing Budgets

| Operation | Target | Notes |
|-----------|--------|-------|
| OP2 read (50 modes) | <2s | pyNastran |
| F06 parse (MEFFMASS) | <0.5s | Regex parsing |
| Modal DB construction | <1s | Data structuring |
| FRF calculation (1000 freq) | <20ms | Vectorized |
| Response PSD | <20ms | Element-wise |
| Spectral integration | <10ms | Trapezoidal |
| Stress recovery (1000 elements) | <200ms | Parallel possible |
| **Total per trial** | **<3s** | Excluding SOL 103 |

### 7.2 Memory Budgets

| Data | Size Estimate | Notes |
|------|---------------|-------|
| Modal DB (100 modes, 10k nodes) | ~50 MB | Full eigenvectors |
| Modal DB (100 modes, 100 nodes) | ~5 MB | Sparse storage |
| Response PSD (1000 freq) | <100 KB | Just arrays |
| Stress shapes (100 modes, 1k elements) | ~50 MB | If needed |

### 7.3 Scalability

Tested configurations:

| Model Size | Modes | Nodes | Elements | Expected Time |
|------------|-------|-------|----------|---------------|
| Small | 50 | 5,000 | 10,000 | <2s |
| Medium | 100 | 50,000 | 100,000 | <5s |
| Large | 200 | 200,000 | 500,000 | <15s |
| Very Large | 500 | 1,000,000 | 2,000,000 | <60s |

---

## 8. Risk Mitigation

### 8.1 Technical Risks

| Risk | Probability | Impact | Mitigation |
|------|-------------|--------|------------|
| pyNastran can't read NX OP2 | Low | High | Test early, have F06 fallback |
| Modal stress not in OP2 | Medium | Medium | Request STRESS(MODAL) output |
| Poor accuracy vs SOL 111 | Low | High | Validate thoroughly, document assumptions |
| Memory issues large models | Medium | Medium | Lazy loading, sparse storage |
| Complex modes (damped) | Low | Low | Detect and warn, use magnitude |

### 8.2 Mitigation Actions

1. **Early OP2 Testing**: Test pyNastran with NX 2412 OP2 files in week 1
2. **F06 Fallback**: Implement F06 parsing for critical data
3. **Validation Suite**: Build comprehensive test cases against SOL 111
4. **Memory Profiling**: Profile with large models, implement streaming if needed
5. **Documentation**: Clear documentation of assumptions and limitations

---

## 9. Future Extensions

### 9.1 Short-Term (3-6 months)

- Sine sweep processor
- Shock response spectrum
- Multi-axis excitation
- Response at multiple nodes
- Automatic mode selection

### 9.2 Medium-Term (6-12 months)

- Fatigue from random vibration
- Acoustic loading (diffuse field)
- Coupled structural-acoustic
- Time-domain transient from PSD
- Design sensitivity (∂response/∂parameter)

### 9.3 Long-Term (12+ months)

- Neural surrogate for modal database
- Real-time response visualization
- Uncertainty quantification
- Multi-physics coupling (thermal + dynamic)
- Foundation model for structural dynamics

---

## 10. Success Metrics

### 10.1 Technical Metrics

- [ ] Grms accuracy <5% vs SOL 111
- [ ] Peak stress accuracy <10% vs SOL 112
- [ ] Processing time <5s per trial
- [ ] Memory <500MB for large models
- [ ] 100% test coverage for core modules

### 10.2 User Metrics

- [ ] Complete optimization in <1 hour (100 trials)
- [ ] Dashboard shows real-time dynamic metrics
- [ ] Documentation enables self-service usage
- [ ] Template covers 80% of use cases

### 10.3 Business Metrics

- [ ] Feature demo ready for marketing video
- [ ] 3+ real-world optimization case studies
- [ ] Positive feedback from beta users
- [ ] Competitive differentiation documented

---

## Appendix A: Reference Equations

### A.1 Modal Superposition

```
Physical response:
u(x,t) = Σᵢ φᵢ(x) · qᵢ(t)

Modal equation of motion:
q̈ᵢ + 2ζᵢωᵢq̇ᵢ + ωᵢ²qᵢ = Lᵢ·a_base(t) / mᵢ

Frequency response function:
Hᵢ(ω) = 1 / (ωᵢ² - ω² + j·2ζᵢωᵢω)
```

### A.2 Random Vibration

```
Response PSD:
Sᵧ(f) = |H(f)|² · Sₓ(f)

Mean square:
σ² = ∫₀^∞ Sᵧ(f) df

Peak (Gaussian):
peak = k · σ,  k ≈ 3 for 3-sigma
```

### A.3 Miles' Equation

```
For SDOF with white noise input:
Grms = √(π/2 · fₙ · Q · PSD)

Where:
- fₙ = natural frequency (Hz)
- Q = 1/(2ζ) = quality factor
- PSD = input PSD level (G²/Hz)
```

---

## Appendix B: File Formats

### B.1 OP2 Data Blocks

| Table | Contains | Used For |
|-------|----------|----------|
| OUGV1 | Eigenvectors | Mode shapes |
| LAMA | Eigenvalues | Frequencies |
| OES1 | Element stress | Modal stress shapes |
| OGPFB1 | Grid point forces | Participation |
| OQMG1 | Modal participation | Effective mass |

### B.2 F06 Tables

| Section | Contains | Parsing Pattern |
|---------|----------|-----------------|
| REAL EIGENVALUES | Frequencies, gen mass | Regex table |
| MODAL EFFECTIVE MASS | Participation | Regex table |
| MODAL PARTICIPATION FACTORS | L factors | Regex table |

### B.3 PSD CSV Format

```csv
frequency_hz,psd_g2_hz
20,0.026
50,0.16
800,0.16
2000,0.026
```

---

*Document End - Dynamic Response Processor Master Plan*