initial commit

.claude/settings.local.json

@@ -0,0 +1,12 @@
+{
+  "permissions": {
+    "allow": [
+      "Bash(if [ -d .cursor/rules ])",
+      "Bash(then ls .cursor/rules)",
+      "Bash(fi)",
+      "Bash(python:*)"
+    ],
+    "deny": [],
+    "ask": []
+  }
+}

CLAUDE.md

@@ -0,0 +1,80 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
+## Project Overview
+
+This is a Python-based transformer modeling toolkit for electrical power transformers. It contains two distinct transformer models and an optimizer:
+
+1. **Simple Ideal Model** ([model.py](model.py)) - Physics-based model for transformers with taps, calculates flux density, losses, and efficiency
+2. **Optimizer** ([optimizer.py](optimizer.py)) - Finds optimal tap settings and input voltage to maximize efficiency for a given load and target power
+3. **Impedance-Based Model** ([old/transformer_model.py](old/transformer_model.py)) - Uses measured open-circuit and short-circuit parameters for AC analysis
+
+## Commands
+
+### Running Simulations
+
+```bash
+# Run the simple transformer model with optimizer example (requires numpy)
+python sim.py
+
+# Run the per-segment resistance example
+python example_per_segment_resistance.py
+
+# Run the impedance-based model example (requires numpy)
+python old/run_model_example.py
+```
+
+### Development Environment
+
+```bash
+# Virtual environment exists at .venv
+# Activate it (Windows):
+.venv\Scripts\activate
+
+# Activate it (Unix/Mac):
+source .venv/bin/activate
+```
+
+## Code Architecture
+
+### Simple Model (model.py)
+
+- **TransformerModel** class: Dataclass-based model with:
+  - Core geometry parameters (Ae_mm2)
+  - Turn counts and tap positions for primary/secondary windings
+  - Copper resistance per turn (supports both uniform and per-segment values)
+  - `simulate()` method: Calculates one operating point given tap number, voltage, frequency, and load
+  - Returns flux density (B_peak_T), currents, powers, and efficiency
+
+**Key concept**: Taps are defined as lists of **incremental turns** for each segment. First element is always 0, subsequent elements specify turns to add. Example: `[0, 128, 23]` means segment 1 adds 128 turns, segment 2 adds 23 turns. Tap numbers are 1-indexed: `tap=1` gives 128 turns, `tap=2` gives 128+23=151 turns total.
+
+**Per-segment resistance**: For windings with different wire gauges per segment, use `primary_Rp_per_turn` and `secondary_Rs_per_turn` lists. Each element specifies R/turn for that segment (e.g., `[0.001, 0.002]` for a 2-segment winding). Legacy single-value `Rp_per_turn`/`Rs_per_turn` still supported.
+
+### Optimizer (optimizer.py)
+
+- **TransformerOptimizer** class: Brute-force search optimizer
+  - `optimize()` method: Searches all tap combinations and voltage sweep to find maximum efficiency
+  - Constraints: Target power delivery (with tolerance), maximum flux density (B_max_T), maximum secondary voltage (Vs_max), input voltage range
+  - Returns **OptimizationResult** with optimal configuration and all operating point details
+  - **Fallback mode** (default): If target power cannot be achieved, finds max achievable power with best efficiency
+
+**Key concept**: For a given load and target power, the optimizer tries all possible tap combinations and input voltages to find the configuration that maximizes efficiency while meeting all constraints (power delivery, flux density, and secondary voltage).
+
+### Impedance Model (old/transformer_model.py)
+
+- **TransformerModel** class: Interpolates measured OC/SC test data
+  - Frequency-dependent Rc and Lm from open-circuit tests
+  - Per-tap series impedance (Req, L_leak) from short-circuit tests
+  - `input_impedance()`: Calculates Zin at frequency/tap/load
+  - `transfer()`: Full AC transfer function with powers and efficiency
+
+**Key concept**: Model uses complex impedance networks with reflected loads through turns ratios.
+
+## Important Implementation Details
+
+- Both models assume **purely resistive loads**
+- Simple model uses **sinusoidal excitation** (4.44*f formula for flux density)
+- Impedance model supports **source resistance** and secondary winding resistance per tap
+- Tap numbering: Both models use 1-indexed taps (tap=1, tap=2, etc.)
+- The `old/` directory contains the earlier impedance-based approach with example sweep functionality

example_per_segment_resistance.py

@@ -0,0 +1,94 @@
+from model import TransformerModel
+from optimizer import TransformerOptimizer
+
+# Example: Transformer with different wire gauges per segment
+# Primary:
+#   - Segment 1 (0 to 50 turns): 22 AWG wire, lower resistance
+#   - Segment 2 (50 to 100 turns): 24 AWG wire, higher resistance
+# Secondary:
+#   - Segment 1 (0 to 150 turns): 28 AWG wire
+#   - Segment 2 (150 to 300 turns): 30 AWG wire, higher resistance
+
+primary_taps = [0, 75, 75]
+secondary_taps = [0, 50, 150, 150, 200]
+
+# Resistance per turn for each segment (example values in ohms/turn)
+# These would be calculated based on wire gauge, length per turn, etc.
+primary_Rp_per_turn = [
+    0.01,
+    0.01,
+]
+
+secondary_Rs_per_turn = [
+    0.004,  
+    0.024,  
+    0.024,  
+    0.024,
+]
+
+tf = TransformerModel(
+    Ae_mm2=354.0,
+    Np_total=150,
+    Ns_total=550,
+    primary_taps=primary_taps,
+    secondary_taps=secondary_taps,
+    primary_Rp_per_turn=primary_Rp_per_turn,
+    secondary_Rs_per_turn=secondary_Rs_per_turn,
+)
+
+
+# Test different tap combinations to see resistance effects
+test_cases = [
+    (1, 1, "Using only first segments (lower R wire)"),
+    (2, 2, "Using both segments (mixed wire gauges)"),
+]
+
+for p_tap, s_tap, description in test_cases:
+    result = tf.simulate(
+        primary_tap=p_tap,
+        secondary_tap=s_tap,
+        Vp_rms=12.0,
+        freq_hz=2000.0,
+        load_ohms=100.0,
+        core_loss_W=0.2,
+    )
+
+    print(f"{description}:")
+    print(f"  Primary tap {p_tap}: {result['Np_eff']} turns")
+    print(f"  Secondary tap {s_tap}: {result['Ns_eff']} turns")
+    print(f"  Output power: {result['P_out_W']:.2f} W")
+    print(f"  Copper loss: {result['P_cu_W']:.2f} W")
+    print(f"  Efficiency: {result['efficiency']*100:.2f}%")
+    print()
+
+# Run optimizer
+print("=" * 70)
+print("OPTIMIZER WITH PER-SEGMENT RESISTANCE")
+print("=" * 70)
+
+opt = TransformerOptimizer(tf)
+
+result_opt = opt.optimize(
+    load_ohms=100.0,
+    target_power_W=10.0,
+    freq_hz=2000.0,
+    Vp_min=5.0,
+    Vp_max=30.0,
+    Vp_step=0.5,
+    B_max_T=0.3,
+    Vs_max=200.0,
+    core_loss_W=0.2,
+    power_tolerance_percent=2.0,
+)
+
+if result_opt:
+    print(f"Optimal configuration:")
+    print(f"  Primary tap: {result_opt.primary_tap} ({result_opt.Np_eff} turns)")
+    print(f"  Secondary tap: {result_opt.secondary_tap} ({result_opt.Ns_eff} turns)")
+    print(f"  Input voltage: {result_opt.Vp_rms:.1f} V")
+    print(f"  Output power: {result_opt.P_out_W:.2f} W")
+    print(f"  Efficiency: {result_opt.efficiency*100:.2f}%")
+    print(f"  Copper loss: {result_opt.P_cu_W:.2f} W")
+    print(f"  Output voltage: {result_opt.Vs_rms:.1f} V")
+else:
+    print("No valid configuration found!")

model.py

@@ -0,0 +1,199 @@
+from dataclasses import dataclass, field
+from typing import List, Tuple, Optional, Dict
+
+
+@dataclass
+class TransformerModel:
+    """
+    Simple transformer model with taps on primary and secondary.
+
+    Assumptions:
+    - Single primary winding with taps.
+    - Single secondary winding with taps.
+    - Ideal magnetic coupling (for now).
+    - Load is purely resistive.
+    """
+    # Core geometry
+    Ae_mm2: float                # effective core area in mm^2
+    # Turn counts (total)
+    Np_total: int
+    Ns_total: int
+    # Tap positions in turns counted from a common reference (e.g. one end of the bobbin)
+    # Must include 0 and total turns if you want to be able to use full winding.
+    primary_taps: List[int] = field(default_factory=lambda: [0])
+    secondary_taps: List[int] = field(default_factory=lambda: [0])
+
+    # Optional: copper resistance per turn (ohms/turn) for each tap segment
+    # If specified, should be a list with length = len(taps) - 1
+    # primary_Rp_per_turn[i] is the R/turn for segment from taps[i] to taps[i+1]
+    # If not specified, falls back to Rp_per_turn/Rs_per_turn
+    primary_Rp_per_turn: Optional[List[float]] = None
+    secondary_Rs_per_turn: Optional[List[float]] = None
+
+    # Legacy: single resistance per turn for all segments (deprecated)
+    Rp_per_turn: float = 0.0
+    Rs_per_turn: float = 0.0
+
+    def _effective_turns(self, taps: List[int], tap: int) -> int:
+        """
+        Compute effective turns for a given tap number.
+
+        tap is 1-indexed: tap=1 uses segment 1, tap=2 uses segments 1+2, etc.
+        Taps list represents incremental turns added at each segment.
+        First element should be 0, subsequent elements are turns to add.
+
+        Example: taps=[0, 50, 150] means:
+        - tap=1: 50 turns (first segment adds 50)
+        - tap=2: 50 + 150 = 200 turns (first segment adds 50, second adds 150)
+        """
+        if tap < 1 or tap >= len(taps):
+            raise ValueError(f"Tap must be between 1 and {len(taps)-1}")
+
+        # Sum all incremental turns: taps[1] + taps[2] + ... + taps[tap]
+        total_turns = sum(taps[1:tap+1])
+
+        return total_turns
+
+    def _resistance_for_tap(
+        self,
+        taps: List[int],
+        tap: int,
+        R_per_turn_list: Optional[List[float]],
+        R_per_turn_legacy: float
+    ) -> float:
+        """
+        Compute total resistance for a given tap.
+
+        Sums resistance across all segments using incremental turns representation.
+        If R_per_turn_list is provided, uses segment-specific resistances.
+        Otherwise uses legacy single R_per_turn value.
+
+        With taps=[0, 50, 150] and tap=2:
+        - Segment 1 has 50 turns
+        - Segment 2 has 150 turns
+        - Total R = R_per_turn[0]*50 + R_per_turn[1]*150
+        """
+        if tap < 1 or tap >= len(taps):
+            raise ValueError(f"Tap must be between 1 and {len(taps)-1}")
+
+        total_resistance = 0.0
+
+        if R_per_turn_list is not None:
+            # Use segment-specific resistance
+            if len(R_per_turn_list) != len(taps) - 1:
+                raise ValueError(f"R_per_turn list must have {len(taps)-1} elements")
+
+            # Sum resistance for each segment (taps[i] represents incremental turns)
+            for i in range(tap):
+                segment_turns = taps[i+1]  # Incremental turns for this segment
+                total_resistance += R_per_turn_list[i] * segment_turns
+        else:
+            # Use legacy single R_per_turn
+            # Total turns is sum of all incremental segments
+            total_turns = sum(taps[1:tap+1])
+            total_resistance = R_per_turn_legacy * total_turns
+
+        return total_resistance
+
+    def simulate(
+        self,
+        primary_tap: int,
+        secondary_tap: int,
+        Vp_rms: float,
+        freq_hz: float,
+        load_ohms: float,
+        core_loss_W: float = 0.0,
+    ) -> Dict[str, float]:
+        """
+        Simulate one operating point.
+
+        Arguments:
+        - primary_tap: tap number (1-indexed), where tap=1 means taps[0] to taps[1], tap=2 means taps[0] to taps[2]
+        - secondary_tap: tap number (1-indexed), where tap=1 means taps[0] to taps[1], tap=2 means taps[0] to taps[2]
+        - Vp_rms: primary RMS voltage (V)
+        - freq_hz: excitation frequency (Hz)
+        - load_ohms: resistive load across the chosen secondary segment (ohms)
+        - core_loss_W: optional core loss at this operating point (W)
+
+        Returns a dict with:
+        - Np_eff, Ns_eff
+        - turns_ratio (Ns_eff / Np_eff)
+        - B_peak_T (Tesla)
+        - Vs_rms, Ip_rms, Is_rms
+        - P_out_W, P_in_W, P_cu_W, P_core_W, efficiency
+        - primary_tap, secondary_tap (echoed back)
+        """
+        if freq_hz <= 0:
+            raise ValueError("Frequency must be > 0")
+        if load_ohms <= 0:
+            raise ValueError("Load resistance must be > 0")
+
+        # Effective turns
+        Np_eff = self._effective_turns(self.primary_taps, primary_tap)
+        Ns_eff = self._effective_turns(self.secondary_taps, secondary_tap)
+
+        if Np_eff == 0 or Ns_eff == 0:
+            raise ValueError("Effective turns cannot be zero")
+
+        # Core area in m^2
+        Ae_m2 = self.Ae_mm2 * 1e-6
+
+        # Peak flux density (sine excitation assumed)
+        B_peak_T = Vp_rms / (4.44 * Np_eff * Ae_m2 * freq_hz)
+
+        # Ideal turn ratio
+        turns_ratio = Ns_eff / Np_eff
+
+        # Ideal secondary RMS voltage
+        Vs_rms_ideal = Vp_rms * turns_ratio
+
+        # Load current and ideal powers (purely resistive)
+        Is_rms_ideal = Vs_rms_ideal / load_ohms
+        Ip_rms_ideal = Is_rms_ideal * turns_ratio  # current ratio inverse of turns
+
+        # Copper resistances for the active segments
+        Rp_seg = self._resistance_for_tap(
+            self.primary_taps,
+            primary_tap,
+            self.primary_Rp_per_turn,
+            self.Rp_per_turn
+        )
+        Rs_seg = self._resistance_for_tap(
+            self.secondary_taps,
+            secondary_tap,
+            self.secondary_Rs_per_turn,
+            self.Rs_per_turn
+        )
+
+        # Copper losses (approx, assuming currents ~ ideal currents)
+        P_cu_p = Ip_rms_ideal**2 * Rp_seg
+        P_cu_s = Is_rms_ideal**2 * Rs_seg
+        P_cu_total = P_cu_p + P_cu_s
+
+        # Output power (resistive load)
+        P_out = Vs_rms_ideal**2 / load_ohms
+
+        # Total input power (approx)
+        P_in = P_out + P_cu_total + core_loss_W
+
+        efficiency = P_out / P_in if P_in > 0 else 1.0
+
+        return {
+            "Np_eff": Np_eff,
+            "Ns_eff": Ns_eff,
+            "turns_ratio": turns_ratio,
+            "B_peak_T": B_peak_T,
+            "Vp_rms": Vp_rms,
+            "Vs_rms": Vs_rms_ideal,
+            "Ip_rms": Ip_rms_ideal,
+            "Is_rms": Is_rms_ideal,
+            "P_out_W": P_out,
+            "P_cu_W": P_cu_total,
+            "P_cu_primary_W": P_cu_p,
+            "P_cu_secondary_W": P_cu_s,
+            "P_core_W": core_loss_W,
+            "P_in_W": P_in,
+            "efficiency": efficiency,
+            "primary_tap": primary_tap,
+            "secondary_tap": secondary_tap,
+        }

optimizer.py

@@ -0,0 +1,240 @@
+from typing import Tuple, Optional, Dict, List
+from dataclasses import dataclass
+from model import TransformerModel
+import itertools
+
+
+@dataclass
+class OptimizationResult:
+    """Result of transformer optimization."""
+    # Optimal configuration
+    primary_tap: int
+    secondary_tap: int
+    Vp_rms: float
+    efficiency: float
+
+    # Operating point details
+    Np_eff: int
+    Ns_eff: int
+    turns_ratio: float
+    B_peak_T: float
+    Vs_rms: float
+    Ip_rms: float
+    Is_rms: float
+    P_out_W: float
+    P_in_W: float
+    P_cu_W: float
+    P_cu_primary_W: float
+    P_cu_secondary_W: float
+    P_core_W: float
+
+    # Constraint satisfaction
+    power_error_percent: float
+    voltage_used: float
+    flux_density_T: float
+
+
+class TransformerOptimizer:
+    """
+    Optimizer for finding best transformer configuration.
+
+    Searches across primary taps, secondary taps, and input voltages to:
+    - Deliver desired power to the load
+    - Maximize efficiency
+    - Stay within flux density and voltage constraints
+    """
+
+    def __init__(self, model: TransformerModel):
+        self.model = model
+
+    def optimize(
+        self,
+        load_ohms: float,
+        target_power_W: float,
+        freq_hz: float,
+        Vp_min: float = 5.0,
+        Vp_max: float = 50.0,
+        Vp_step: float = 0.5,
+        B_max_T: float = 0.3,
+        Vs_max: float = float('inf'),
+        core_loss_W: float = 0.0,
+        power_tolerance_percent: float = 2.0,
+        fallback_max_power: bool = True,
+    ) -> Optional[OptimizationResult]:
+        """
+        Find optimal transformer configuration for given constraints.
+
+        Arguments:
+        - load_ohms: Resistive load (ohms)
+        - target_power_W: Desired output power (W)
+        - freq_hz: Operating frequency (Hz)
+        - Vp_min, Vp_max, Vp_step: Input voltage search range
+        - B_max_T: Maximum allowed flux density (Tesla)
+        - Vs_max: Maximum allowed secondary voltage (V)
+        - core_loss_W: Core loss at this operating point
+        - power_tolerance_percent: Acceptable power delivery error (%)
+        - fallback_max_power: If True and target cannot be met, find max power with best efficiency
+
+        Returns:
+        - OptimizationResult with best configuration, or None if no valid solution
+        """
+
+        if load_ohms <= 0:
+            raise ValueError("Load must be > 0")
+        if target_power_W <= 0:
+            raise ValueError("Target power must be > 0")
+        if Vp_min >= Vp_max:
+            raise ValueError("Vp_min must be < Vp_max")
+
+        # Generate all possible tap numbers
+        primary_taps = list(range(1, len(self.model.primary_taps)))
+        secondary_taps = list(range(1, len(self.model.secondary_taps)))
+
+        if not primary_taps or not secondary_taps:
+            raise ValueError("Need at least 2 taps on each winding")
+
+        best_result = None
+        best_efficiency = -1.0
+
+        # Fallback tracking: find max achievable power under target
+        fallback_result = None
+        fallback_max_power_achieved = -1.0
+        fallback_best_efficiency = -1.0
+
+        # Search over all tap combinations
+        for p_tap, s_tap in itertools.product(primary_taps, secondary_taps):
+            # For each tap combination, find the voltage that delivers target power
+            result = self._optimize_voltage_for_tap(
+                primary_tap=p_tap,
+                secondary_tap=s_tap,
+                load_ohms=load_ohms,
+                target_power_W=target_power_W,
+                freq_hz=freq_hz,
+                Vp_min=Vp_min,
+                Vp_max=Vp_max,
+                Vp_step=Vp_step,
+                B_max_T=B_max_T,
+                Vs_max=Vs_max,
+                core_loss_W=core_loss_W,
+                power_tolerance_percent=power_tolerance_percent,
+                fallback_max_power=fallback_max_power,
+            )
+
+            if result is not None:
+                power_error = abs(result.P_out_W - target_power_W) / target_power_W * 100
+
+                # Check if this meets the target power tolerance
+                if power_error <= power_tolerance_percent:
+                    if result.efficiency > best_efficiency:
+                        best_efficiency = result.efficiency
+                        best_result = result
+                # Track fallback: maximize power below target, then maximize efficiency
+                elif fallback_max_power and result.P_out_W < target_power_W:
+                    if (result.P_out_W > fallback_max_power_achieved or
+                        (result.P_out_W == fallback_max_power_achieved and result.efficiency > fallback_best_efficiency)):
+                        fallback_max_power_achieved = result.P_out_W
+                        fallback_best_efficiency = result.efficiency
+                        fallback_result = result
+
+        # Return best result if found, otherwise fallback
+        return best_result if best_result is not None else fallback_result
+
+    def _optimize_voltage_for_tap(
+        self,
+        primary_tap: int,
+        secondary_tap: int,
+        load_ohms: float,
+        target_power_W: float,
+        freq_hz: float,
+        Vp_min: float,
+        Vp_max: float,
+        Vp_step: float,
+        B_max_T: float,
+        Vs_max: float,
+        core_loss_W: float,
+        power_tolerance_percent: float,
+        fallback_max_power: bool,
+    ) -> Optional[OptimizationResult]:
+        """
+        For a given tap configuration, find the best voltage.
+
+        In normal mode: returns highest efficiency that meets target power ± tolerance
+        In fallback mode: returns highest efficiency at max achievable power below target
+
+        Returns OptimizationResult if valid solution found, else None.
+        """
+        best_within_tolerance = None
+        best_efficiency_within_tolerance = -1.0
+
+        # For fallback: track max power achievable and best efficiency at that power
+        max_power_below_target = -1.0
+        best_at_max_power = None
+        best_efficiency_at_max_power = -1.0
+
+        # Generate voltage sweep
+        import numpy as np
+        voltages = np.arange(Vp_min, Vp_max + Vp_step/2, Vp_step)
+
+        for Vp_rms in voltages:
+            try:
+                sim = self.model.simulate(
+                    primary_tap=primary_tap,
+                    secondary_tap=secondary_tap,
+                    Vp_rms=float(Vp_rms),
+                    freq_hz=freq_hz,
+                    load_ohms=load_ohms,
+                    core_loss_W=core_loss_W,
+                )
+            except (ValueError, ZeroDivisionError):
+                continue
+
+            # Check constraints
+            if sim["B_peak_T"] > B_max_T:
+                continue
+
+            if sim["Vs_rms"] > Vs_max:
+                continue
+
+            power_error = abs(sim["P_out_W"] - target_power_W) / target_power_W * 100
+
+            result = OptimizationResult(
+                primary_tap=primary_tap,
+                secondary_tap=secondary_tap,
+                Vp_rms=float(Vp_rms),
+                efficiency=sim["efficiency"],
+                Np_eff=sim["Np_eff"],
+                Ns_eff=sim["Ns_eff"],
+                turns_ratio=sim["turns_ratio"],
+                B_peak_T=sim["B_peak_T"],
+                Vs_rms=sim["Vs_rms"],
+                Ip_rms=sim["Ip_rms"],
+                Is_rms=sim["Is_rms"],
+                P_out_W=sim["P_out_W"],
+                P_in_W=sim["P_in_W"],
+                P_cu_W=sim["P_cu_W"],
+                P_cu_primary_W=sim["P_cu_primary_W"],
+                P_cu_secondary_W=sim["P_cu_secondary_W"],
+                P_core_W=sim["P_core_W"],
+                power_error_percent=power_error,
+                voltage_used=float(Vp_rms),
+                flux_density_T=sim["B_peak_T"],
+            )
+
+            # Track solutions within tolerance
+            if power_error <= power_tolerance_percent:
+                if sim["efficiency"] > best_efficiency_within_tolerance:
+                    best_efficiency_within_tolerance = sim["efficiency"]
+                    best_within_tolerance = result
+
+            # Track fallback: max power below target with best efficiency
+            if fallback_max_power and sim["P_out_W"] < target_power_W:
+                if sim["P_out_W"] > max_power_below_target:
+                    max_power_below_target = sim["P_out_W"]
+                    best_efficiency_at_max_power = sim["efficiency"]
+                    best_at_max_power = result
+                elif sim["P_out_W"] == max_power_below_target and sim["efficiency"] > best_efficiency_at_max_power:
+                    best_efficiency_at_max_power = sim["efficiency"]
+                    best_at_max_power = result
+
+        # Return best within tolerance if found, otherwise fallback
+        return best_within_tolerance if best_within_tolerance is not None else best_at_max_power

sim.py

@@ -0,0 +1,121 @@
+from model import TransformerModel
+from optimizer import TransformerOptimizer
+
+# Example: E55/28/21 N27 transformer with one primary tap and simple secondary
+
+primary_taps = [0, 75, 75]
+secondary_taps = [0, 100, 150, 150, 150]
+
+# Resistance per turn for each segment (example values in ohms/turn)
+# These would be calculated based on wire gauge, length per turn, etc.
+primary_Rp_per_turn = [
+    0.01,
+    0.01,
+]
+
+secondary_Rs_per_turn = [
+    0.004,  
+    0.024,  
+    0.024,  
+    0.024,
+]
+
+tf = TransformerModel(
+    Ae_mm2=354.0,
+    Np_total=150,
+    Ns_total=550,
+    primary_taps=primary_taps,
+    secondary_taps=secondary_taps,
+    primary_Rp_per_turn=primary_Rp_per_turn,
+    secondary_Rs_per_turn=secondary_Rs_per_turn,
+)
+
+
+# (Optional) set copper resistance per turn when you know it
+# For example, from earlier rough calc:
+#   Rp_total ~ 0.9 Ω for 151 turns -> ~0.9/151 per turn
+#   Rs_total ~ 10.8 Ω for 567 turns -> ~10.8/567 per turn
+# tf.Rp_per_turn = 0.9 / 151
+# tf.Rs_per_turn = 10.8 / 567
+
+
+print("=" * 70)
+print("MANUAL SIMULATION EXAMPLE")
+print("=" * 70)
+
+# Use tap 1 for primary (taps[0] to taps[1] => 0 to 128 turns)
+# Use tap 1 for secondary (taps[0] to taps[1] => 0 to 567 turns)
+primary_tap = 2
+secondary_tap = 1
+
+# Say you apply 12 Vrms at 2 kHz into the primary,
+# and you choose a load such that power is drawn from the secondary.
+result_manual = tf.simulate(
+    primary_tap=primary_tap,
+    secondary_tap=secondary_tap,
+    Vp_rms=18,
+    freq_hz=256.0,
+    load_ohms=6.0,
+    core_loss_W=0.3,  # example; refine later with measurements
+)
+
+print(f"Manual configuration:")
+print(f"  Primary tap: {primary_tap} ({result_manual['Np_eff']} turns)")
+print(f"  Secondary tap: {secondary_tap} ({result_manual['Ns_eff']} turns)")
+print(f"  Input voltage: {result_manual['Vp_rms']:.1f} V")
+print(f"  Turns ratio: {result_manual['turns_ratio']:.2f}")
+print(f"  Output power: {result_manual['P_out_W']:.2f} W")
+print(f"  Efficiency: {result_manual['efficiency']*100:.2f}%")
+print(f"  Flux density: {result_manual['B_peak_T']:.3f} T")
+print(f"  Input current: {result_manual['Ip_rms']:.3f} A")
+print(f"  Output voltage: {result_manual['Vs_rms']:.1f} V")
+print(f"  Output current: {result_manual['Is_rms']:.3f} A")
+print(f"  Copper loss (total): {result_manual['P_cu_W']:.2f} W")
+print(f"    Primary: {result_manual['P_cu_primary_W']:.2f} W")
+print(f"    Secondary: {result_manual['P_cu_secondary_W']:.2f} W")
+print(f"  Core loss: {result_manual['P_core_W']:.2f} W")
+print()
+
+
+print("=" * 70)
+print("OPTIMIZER EXAMPLE")
+print("=" * 70)
+
+# Create optimizer
+opt = TransformerOptimizer(tf)
+
+# Find optimal configuration for delivering 10W to a 1000 ohm load at 2 kHz
+result_opt = opt.optimize(
+    load_ohms=1500.0,
+    target_power_W=10.0,
+    freq_hz=45000.0,
+    Vp_min=1.0,
+    Vp_max=36.0,
+    Vp_step=0.5,
+    B_max_T=0.4,  # maximum flux density constraint
+    Vs_max=130.0,  # maximum secondary voltage constraint
+    core_loss_W=0.3,
+    power_tolerance_percent=2.0,  # accept ±2% power delivery error
+)
+
+if result_opt:
+    target_power = 10.0
+    if result_opt.power_error_percent > 2.0:
+        print(f"NOTE: Target {target_power}W not achievable. Showing max power configuration.")
+    print(f"Optimal configuration found:")
+    print(f"  Primary tap: {result_opt.primary_tap} ({result_opt.Np_eff} turns)")
+    print(f"  Secondary tap: {result_opt.secondary_tap} ({result_opt.Ns_eff} turns)")
+    print(f"  Input voltage: {result_opt.Vp_rms:.1f} V")
+    print(f"  Turns ratio: {result_opt.turns_ratio:.2f}")
+    print(f"  Output power: {result_opt.P_out_W:.2f} W (error: {result_opt.power_error_percent:.2f}%)")
+    print(f"  Efficiency: {result_opt.efficiency*100:.2f}%")
+    print(f"  Flux density: {result_opt.B_peak_T:.3f} T")
+    print(f"  Input current: {result_opt.Ip_rms:.3f} A")
+    print(f"  Output voltage: {result_opt.Vs_rms:.1f} V")
+    print(f"  Output current: {result_opt.Is_rms:.3f} A")
+    print(f"  Copper loss (total): {result_opt.P_cu_W:.2f} W")
+    print(f"    Primary: {result_opt.P_cu_primary_W:.2f} W")
+    print(f"    Secondary: {result_opt.P_cu_secondary_W:.2f} W")
+    print(f"  Core loss: {result_opt.P_core_W:.2f} W")
+else:
+    print("No valid configuration found within constraints!")