adding tordoidal simulator and drawing

README_WEBAPP.md

@@ -0,0 +1,83 @@
+# Transformer Optimizer Web App
+
+A Flask-based web application for optimizing transformer designs and running manual simulations.
+
+## Pages
+
+### 1. Optimizer (http://localhost:5000/)
+Automatically finds the best tap configuration and input voltage to maximize efficiency.
+
+### 2. Manual Simulation (http://localhost:5000/manual)
+Directly specify tap settings and input voltage to see performance results.
+
+## Features
+
+### Optimizer Page
+- **Interactive Sliders** for all parameters:
+  - Load resistance (5Ω to 10kΩ)
+  - Target power output (1W to 100W)
+  - Operating frequency (100Hz to 50kHz)
+  - Core loss (0W to 5W)
+  - Maximum flux density (0.05T to 0.5T)
+  - Maximum secondary voltage (10V to 500V)
+  - Input voltage range (min/max)
+
+- **Real-time Optimization** using the transformer optimizer
+
+### Manual Simulation Page
+- **Tap Selection**: Choose primary and secondary taps from dropdowns
+- **Input Parameters**: Set input voltage and frequency with sliders
+- **Load Conditions**: Configure load resistance and core loss
+
+### Both Pages Display:
+- Efficiency (color-coded: green >90%, yellow >80%, red <80%)
+- Tap selections with turn counts
+- All voltages and currents
+- Separate copper losses for primary and secondary
+- Flux density and turns ratio
+- Input and output power
+
+## Installation
+
+1. Install dependencies:
+```bash
+pip install -r requirements.txt
+```
+
+## Running the App
+
+```bash
+python app.py
+```
+
+Then open your browser to: http://localhost:5000
+
+## Configuration
+
+Edit `app.py` and modify the `get_default_transformer()` function to match your transformer specifications:
+- Tap positions
+- Core area
+- Turn counts
+- Resistance per turn for each segment
+
+## API Endpoints
+
+- `GET /` - Optimizer page
+- `GET /manual` - Manual simulation page
+- `POST /api/optimize` - Run optimization with parameters
+- `POST /api/simulate` - Run manual simulation
+- `GET /api/transformer_info` - Get transformer configuration
+
+## Usage
+
+### Optimizer Mode
+1. Adjust sliders to set your desired operating conditions
+2. Click "⚡ Optimize Transformer"
+3. View the optimized configuration and performance metrics
+4. The optimizer will find the best tap combination and input voltage to maximize efficiency while meeting your constraints
+
+### Manual Simulation Mode
+1. Select primary and secondary taps from dropdowns
+2. Adjust input voltage, frequency, load, and core loss with sliders
+3. Click "▶️ Run Simulation"
+4. View the performance results for your specified configuration

app.py

@@ -0,0 +1,197 @@
+from flask import Flask, render_template, request, jsonify
+from model import TransformerModel
+from optimizer import TransformerOptimizer
+
+app = Flask(__name__)
+
+# Default transformer configuration
+# You can modify this to match your actual transformer
+def get_default_transformer():
+
+    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,
+        Ve_mm3=43900.0,
+        use_core_loss_model=True,
+        Np_total=150,
+        Ns_total=250,
+        primary_taps=primary_taps,
+        secondary_taps=secondary_taps,
+        primary_Rp_per_turn=primary_Rp_per_turn,
+        secondary_Rs_per_turn=secondary_Rs_per_turn,
+    )
+    return tf
+
+
+@app.route('/')
+def index():
+    return render_template('index.html')
+
+
+@app.route('/manual')
+def manual():
+    return render_template('manual.html')
+
+
+@app.route('/api/simulate', methods=['POST'])
+def simulate():
+    """Run manual simulation with specified tap and voltage"""
+    try:
+        data = request.json
+
+        # Extract parameters
+        primary_tap = int(data.get('primary_tap', 1))
+        secondary_tap = int(data.get('secondary_tap', 1))
+        Vp_rms = float(data.get('Vp_rms', 12))
+        freq_hz = float(data.get('freq_hz', 2000))
+        load_ohms = float(data.get('load_ohms', 100))
+
+        # Create transformer
+        tf = get_default_transformer()
+
+        # Run simulation (core loss is calculated automatically)
+        result = tf.simulate(
+            primary_tap=primary_tap,
+            secondary_tap=secondary_tap,
+            Vp_rms=Vp_rms,
+            freq_hz=freq_hz,
+            load_ohms=load_ohms,
+            core_loss_W=0.0,  # Will be calculated by model
+        )
+
+        return jsonify({
+            'success': True,
+            'result': {
+                'primary_tap': result['primary_tap'],
+                'secondary_tap': result['secondary_tap'],
+                'Np_eff': result['Np_eff'],
+                'Ns_eff': result['Ns_eff'],
+                'Vp_rms': round(result['Vp_rms'], 2),
+                'Vs_rms': round(result['Vs_rms'], 2),
+                'Ip_rms': round(result['Ip_rms'], 3),
+                'Is_rms': round(result['Is_rms'], 3),
+                'turns_ratio': round(result['turns_ratio'], 3),
+                'P_out_W': round(result['P_out_W'], 2),
+                'P_in_W': round(result['P_in_W'], 2),
+                'P_cu_W': round(result['P_cu_W'], 2),
+                'P_cu_primary_W': round(result['P_cu_primary_W'], 3),
+                'P_cu_secondary_W': round(result['P_cu_secondary_W'], 3),
+                'P_core_W': round(result['P_core_W'], 2),
+                'efficiency': round(result['efficiency'] * 100, 2),
+                'B_peak_T': round(result['B_peak_T'], 4),
+            }
+        })
+
+    except Exception as e:
+        return jsonify({
+            'success': False,
+            'error': str(e)
+        }), 400
+
+
+@app.route('/api/optimize', methods=['POST'])
+def optimize():
+    try:
+        data = request.json
+
+        # Extract parameters
+        load_ohms = float(data.get('load_ohms', 100))
+        target_power_W = float(data.get('target_power_W', 10))
+        freq_hz = float(data.get('freq_hz', 2000))
+        Vp_min = float(data.get('Vp_min', 5))
+        Vp_max = float(data.get('Vp_max', 36))
+        Vp_step = float(data.get('Vp_step', 0.5))
+        B_max_T = float(data.get('B_max_T', 0.3))
+        Vs_max = float(data.get('Vs_max', 200))
+        Is_max = float(data.get('Is_max', 1.5))
+        power_tolerance_percent = float(data.get('power_tolerance_percent', 2.0))
+
+        # Create transformer and optimizer
+        tf = get_default_transformer()
+        opt = TransformerOptimizer(tf)
+
+        # Run optimization (core loss is calculated automatically)
+        result = opt.optimize(
+            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,
+            Is_max=Is_max,
+            core_loss_W=0.0,  # Will be calculated by model
+            power_tolerance_percent=power_tolerance_percent,
+        )
+
+        if result:
+            return jsonify({
+                'success': True,
+                'result': {
+                    'primary_tap': result.primary_tap,
+                    'secondary_tap': result.secondary_tap,
+                    'Np_eff': result.Np_eff,
+                    'Ns_eff': result.Ns_eff,
+                    'Vp_rms': round(result.Vp_rms, 2),
+                    'Vs_rms': round(result.Vs_rms, 2),
+                    'Ip_rms': round(result.Ip_rms, 3),
+                    'Is_rms': round(result.Is_rms, 3),
+                    'turns_ratio': round(result.turns_ratio, 3),
+                    'P_out_W': round(result.P_out_W, 2),
+                    'P_in_W': round(result.P_in_W, 2),
+                    'P_cu_W': round(result.P_cu_W, 2),
+                    'P_cu_primary_W': round(result.P_cu_primary_W, 3),
+                    'P_cu_secondary_W': round(result.P_cu_secondary_W, 3),
+                    'P_core_W': round(result.P_core_W, 2),
+                    'efficiency': round(result.efficiency * 100, 2),
+                    'B_peak_T': round(result.B_peak_T, 4),
+                    'power_error_percent': round(result.power_error_percent, 2),
+                }
+            })
+        else:
+            return jsonify({
+                'success': False,
+                'error': 'No valid configuration found within constraints'
+            })
+
+    except Exception as e:
+        return jsonify({
+            'success': False,
+            'error': str(e)
+        }), 400
+
+
+@app.route('/api/transformer_info', methods=['GET'])
+def transformer_info():
+    """Return transformer configuration information"""
+    tf = get_default_transformer()
+    return jsonify({
+        'primary_taps': tf.primary_taps,
+        'secondary_taps': tf.secondary_taps,
+        'Ae_mm2': tf.Ae_mm2,
+        'Np_total': tf.Np_total,
+        'Ns_total': tf.Ns_total,
+    })
+
+
+if __name__ == '__main__':
+    app.run(debug=True, host='0.0.0.0', port=5000)

core.py

@@ -0,0 +1,207 @@
+import numpy as np
+
+# -------------------------------------------------------------------
+# 1. Put your digitized data here
+#    Units:
+#      - B in Tesla
+#      - frequency in kHz (x-axis of your plot)
+#      - Pv in kW/m^3 (y-axis of your plot)
+#
+# IMPORTANT:
+#   The numbers below are ONLY EXAMPLES / PLACEHOLDERS.
+#   Replace each "f_khz" and "Pv_kW_m3" array with points
+#   read from your datasheet plot (solid 25 °C lines).
+# -------------------------------------------------------------------
+
+core_loss_data_25C = {
+    0.025: {  # 25 mT
+        "f_khz": np.array([2.93084374470676,
+            8.7361542468832,
+            44.706096185109,
+            125.963436020903,
+            193.225112880907,
+            ], dtype=float),
+        "Pv_kW_m3": np.array([0.177193262696068,
+            0.679670055798598,
+            5.03342212801684,
+            17.4677120801992,
+            30.0788251804309,
+            ], dtype=float),
+    },
+    0.050: {  # 50 mT
+        "f_khz": np.array([2.9975926563905,
+            9.24208952674728,
+            38.618648046052,
+            188.92247157063,
+            ], dtype=float),
+        "Pv_kW_m3": np.array([0.917768017464649,
+            3.47034515205292,
+            19.0328184767803,
+            123.927503749051,
+            ], dtype=float),
+    },
+    0.100: {  # 100 mT
+        "f_khz": np.array([2.48,
+            8.84,
+            34.12,
+            87.85,
+            193.23,
+            ], dtype=float),
+        "Pv_kW_m3": np.array([3.14,
+            13.89,
+            69.94,
+            216.47,
+            564.36,
+            ], dtype=float),
+    },
+    0.200: {  # 200 mT
+        "f_khz": np.array([2.14,
+                8.64,
+                32.25,
+                83.05,
+                195.41,
+                ], dtype=float),
+        "Pv_kW_m3": np.array([10.44,
+            56.44,
+            276.05,
+            891.89,
+            2497.56,
+            ], dtype=float),
+    },
+    0.300: {  # 300 mT
+        "f_khz": np.array([2.19,
+            5.03,
+            12.39,
+            24.89,
+            49.47,
+            ], dtype=float),
+        "Pv_kW_m3": np.array([26.83,
+            67.97,
+            193.07,
+            442.55,
+            1000,
+            ], dtype=float),
+    },
+}
+
+
+# -------------------------------------------------------------------
+# 2. Helper: 1-D log-log interpolation in frequency
+# -------------------------------------------------------------------
+
+def _loglog_interp(x, xp, fp):
+    """
+    Log-log interpolation: fp vs xp, evaluated at x.
+
+    All arguments are assumed > 0.
+    """
+    x = np.asarray(x, dtype=float)
+    xp = np.asarray(xp, dtype=float)
+    fp = np.asarray(fp, dtype=float)
+
+    logx = np.log10(x)
+    logxp = np.log10(xp)
+    logfp = np.log10(fp)
+
+    logy = np.interp(logx, logxp, logfp)
+    return 10.0 ** logy
+
+
+# -------------------------------------------------------------------
+# 3. Main function: Pv(B, f) interpolation
+# -------------------------------------------------------------------
+
+def core_loss_Pv(B_T, f_Hz, data=core_loss_data_25C):
+    """
+    Interpolate core loss density P_v for given flux density and frequency.
+
+    Parameters
+    ----------
+    B_T : float
+        Peak flux density in Tesla (e.g. 0.05 for 50 mT).
+    f_Hz : float or array_like
+        Frequency in Hz.
+    data : dict
+        Dict mapping B_T -> {"f_khz": array, "Pv_kW_m3": array}.
+
+    Returns
+    -------
+    Pv_kW_m3 : float or ndarray
+        Core loss density in kW/m^3 at 25 °C.
+    """
+    B_T = float(B_T)
+    f_khz = np.asarray(f_Hz, dtype=float) / 1e3  # plot uses kHz
+
+    if np.any(f_khz <= 0):
+        raise ValueError("Frequency must be > 0 Hz")
+
+    # Sorted list of available B values from your digitized curves
+    B_list = np.array(sorted(data.keys()))
+
+    # Handle out-of-range B values
+    if B_T < B_list.min():
+        raise ValueError(
+            f"B={B_T:.3f} T is below available range "
+            f"[{B_list.min():.3f}, {B_list.max():.3f}] T"
+        )
+
+    # For B > max data (0.3T), extrapolate up to 0.4T using two highest curves
+    # If B > 0.4T, clamp to the extrapolated value at 0.4T
+    if B_T > B_list.max():
+        # Use the two highest B values for extrapolation
+        B1 = B_list[-2]  # 0.2T
+        B2 = B_list[-1]  # 0.3T
+
+        d1 = data[B1]
+        d2 = data[B2]
+
+        # Interpolate along frequency (log-log) on each B curve
+        Pv1 = _loglog_interp(f_khz, d1["f_khz"], d1["Pv_kW_m3"])
+        Pv2 = _loglog_interp(f_khz, d2["f_khz"], d2["Pv_kW_m3"])
+
+        # Extrapolate in log(Pv) vs B (Pv ~ B^b approximately)
+        logPv1 = np.log10(Pv1)
+        logPv2 = np.log10(Pv2)
+
+        # Clamp B_T to 0.4T for extrapolation calculation
+        B_T_clamped = min(B_T, 0.4)
+        alpha = (B_T_clamped - B1) / (B2 - B1)
+        logPv = (1.0 - alpha) * logPv1 + alpha * logPv2
+
+        return 10.0 ** logPv
+
+    # If B_T is (approximately) one of the tabulated curves, just interpolate in f
+    idx_match = np.where(np.isclose(B_list, B_T, rtol=1e-4, atol=1e-6))[0]
+    if idx_match.size:
+        B_key = B_list[idx_match[0]]
+        d = data[B_key]
+        return _loglog_interp(f_khz, d["f_khz"], d["Pv_kW_m3"])
+
+    # Otherwise, find the two B curves that bracket B_T
+    idx = np.searchsorted(B_list, B_T)
+    B1 = B_list[idx - 1]
+    B2 = B_list[idx]
+
+    d1 = data[B1]
+    d2 = data[B2]
+
+    # Interpolate along frequency (log-log) on each B curve
+    Pv1 = _loglog_interp(f_khz, d1["f_khz"], d1["Pv_kW_m3"])
+    Pv2 = _loglog_interp(f_khz, d2["f_khz"], d2["Pv_kW_m3"])
+
+    # Now interpolate between B1 and B2.
+    # Do it in log(Pv) vs B (Pv ~ B^b approximately).
+    logPv1 = np.log10(Pv1)
+    logPv2 = np.log10(Pv2)
+
+    alpha = (B_T - B1) / (B2 - B1)
+    logPv = (1.0 - alpha) * logPv1 + alpha * logPv2
+
+    return 10.0 ** logPv
+
+
+
+if __name__ == '__main__':
+    p = core_loss_Pv(0.2, 25_000)
+    print (p)
+

designer.py

@@ -0,0 +1,641 @@
+"""
+Toroid Winding Designer
+=======================
+Determines the feasibility of winding a toroidal transformer given core
+geometry, wire gauge, and turn counts.  For each winding segment (tap) it
+reports:
+
+  - Whether the turns physically fit on the core
+  - Number of complete layers required and turns per layer
+  - Total wire length
+  - DC resistance
+  - Wire weight
+
+Wire gauge table
+----------------
+Diameters in the AWG table already include insulation (heavy-build magnet
+wire), so no additional insulation allowance is applied.
+
+Geometry assumptions (toroid cross-section is a rectangle w × h)
+-----------------------------------------------------------------
+Winding builds simultaneously on ALL four faces of the rectangular
+cross-section:
+  - Inner bore radius shrinks  by d_wire per layer
+  - Outer radius grows         by d_wire per layer
+  - Available winding height   shrinks by 2*d_wire per layer
+
+For layer n (0-indexed) the per-turn path length is:
+  L_turn(n) = π*(OD/2 + n*d + ID/2 - (n+1)*d) + 2*(h - 2*n*d)
+            = π*((OD + ID)/2 - d) + 2*h - 4*n*d
+
+Turns that fit in layer n (minimum of inner-bore and height constraints):
+  turns_inner(n) = floor( π*(ID - 2*(n+1)*d) / d )
+  turns_height(n) = floor( (h - 2*n*d) / d )
+  turns_per_layer(n) = min(turns_inner(n), turns_height(n))
+
+Winding stops when either constraint hits zero (no more room).
+
+The fill_factor parameter (default 0.35) limits the fraction of the
+available window area that can be used.  Window area = π*(ID/2)² for a
+toroid.  Wire cross-section area per turn = π*(d/2)².  Max turns from
+fill factor = floor(fill_factor * window_area / wire_area).
+"""
+
+from __future__ import annotations
+
+import math
+from dataclasses import dataclass, field
+from typing import Optional
+
+# ---------------------------------------------------------------------------
+# Wire gauge table
+# AWG: (diameter_mm including insulation, area_mm2 conductor, ohm_per_km)
+# ---------------------------------------------------------------------------
+
+AWG_TABLE: dict[int, tuple[float, float, float]] = {
+    20: (0.813, 0.518, 33.292),
+    21: (0.724, 0.410, 41.984),
+    22: (0.645, 0.326, 52.939),
+    23: (0.574, 0.258, 66.781),
+    24: (0.511, 0.205, 84.198),
+    25: (0.455, 0.162, 106.17),
+    26: (0.404, 0.129, 133.86),
+    27: (0.361, 0.102, 168.62),
+    28: (0.320, 0.081, 212.87),
+    29: (0.287, 0.064, 268.4),
+    30: (0.254, 0.051, 338.5),
+}
+
+# Copper density in g/mm³
+_COPPER_DENSITY_G_MM3 = 8.96e-3
+
+# ---------------------------------------------------------------------------
+# Wire cost table
+# AWG: (spool_length_m, cost_usd_per_spool, part_number)
+# ---------------------------------------------------------------------------
+
+WIRE_COST_TABLE: dict[int, tuple[float, float, str]] = {
+    20: (972,  230.00, "20SNS10"),
+    22: (1545, 238.23, "22SNS10"),
+    24: (2447, 222.53, "24SNS10"),
+    28: (6178, 224.00, "28SNS10"),
+}
+
+
+@dataclass
+class WireSpec:
+    """Properties of a single AWG wire."""
+    awg: int
+    diameter_mm: float    # including insulation
+    area_mm2: float       # conductor cross-section
+    ohm_per_km: float     # DC resistance
+    spool_length_m: Optional[float] = None   # None if not in cost table
+    cost_per_spool: Optional[float] = None   # USD
+    part_number: Optional[str] = None
+
+    @property
+    def radius_mm(self) -> float:
+        return self.diameter_mm / 2.0
+
+    @property
+    def cost_per_m(self) -> Optional[float]:
+        """USD per metre, or None if cost data unavailable."""
+        if self.spool_length_m and self.cost_per_spool:
+            return self.cost_per_spool / self.spool_length_m
+        return None
+
+    @classmethod
+    def from_awg(cls, awg: int) -> "WireSpec":
+        if awg not in AWG_TABLE:
+            raise ValueError(f"AWG {awg} not in table. Available: {sorted(AWG_TABLE)}")
+        d, a, r = AWG_TABLE[awg]
+        spool_m = cost_usd = pn = None
+        if awg in WIRE_COST_TABLE:
+            spool_m, cost_usd, pn = WIRE_COST_TABLE[awg]
+        return cls(awg=awg, diameter_mm=d, area_mm2=a, ohm_per_km=r,
+                   spool_length_m=spool_m, cost_per_spool=cost_usd, part_number=pn)
+
+
+@dataclass
+class ToroidCore:
+    """
+    Toroidal core geometry and magnetic material parameters.
+
+    Parameters
+    ----------
+    ID_mm : float
+        Inner diameter of the core in mm.
+    OD_mm : float
+        Outer diameter of the core in mm.
+    height_mm : float
+        Height (axial length) of the core in mm.
+    Ae_mm2 : float, optional
+        Effective magnetic cross-section area (mm²).  Used for flux-density
+        calculation.  Defaults to cross_section_area_mm2 (geometric value).
+    Ve_mm3 : float, optional
+        Effective core volume (mm³).  Used for core loss calculation.
+        Defaults to Ae_mm2 * mean_path_length_mm.  Set to 0 to disable
+        core loss.
+    Pv_func : callable, optional
+        Core loss density function with signature::
+
+            Pv_func(f_hz: float, B_T: float) -> float   # kW / m³
+
+        If None and Ve_mm3 > 0, the simulator falls back to the built-in
+        ``core.core_loss_Pv`` interpolation (which also returns kW/m³).
+        Pass ``lambda f, B: 0.0`` to disable core loss while keeping a
+        non-zero volume.
+    """
+    ID_mm: float
+    OD_mm: float
+    height_mm: float
+    Ae_mm2: Optional[float] = None
+    Ve_mm3: Optional[float] = None
+    Pv_func: Optional[object] = None   # Callable[[float, float], float] | None
+
+    def __post_init__(self):
+        if self.ID_mm <= 0 or self.OD_mm <= self.ID_mm or self.height_mm <= 0:
+            raise ValueError("Require 0 < ID < OD and height > 0")
+        if self.Ae_mm2 is not None and self.Ae_mm2 <= 0:
+            raise ValueError("Ae_mm2 must be > 0")
+        if self.Ve_mm3 is not None and self.Ve_mm3 < 0:
+            raise ValueError("Ve_mm3 must be >= 0")
+
+    @property
+    def window_area_mm2(self) -> float:
+        """Area of the inner bore (the winding window)."""
+        return math.pi * (self.ID_mm / 2.0) ** 2
+
+    @property
+    def cross_section_area_mm2(self) -> float:
+        """Rectangular cross-section area of the core material."""
+        return ((self.OD_mm - self.ID_mm) / 2.0) * self.height_mm
+
+    @property
+    def mean_path_length_mm(self) -> float:
+        """Mean magnetic path length through the core centre."""
+        return math.pi * (self.OD_mm + self.ID_mm) / 2.0
+
+
+@dataclass
+class WindingSpec:
+    """
+    Specification for one winding (primary or secondary).
+
+    Parameters
+    ----------
+    awg : int or list[int]
+        Wire gauge(s).  A single int applies to all segments.  A list must
+        have one entry per segment (i.e. len(taps) - 1 entries).
+        Example: awg=26  or  awg=[26, 28]  for a 2-segment winding.
+    taps : list[int]
+        Incremental turns per segment, same convention as model.py.
+        First element is always 0; each subsequent element is the additional
+        turns added by that tap segment.
+        Example: [0, 128, 23]  → tap 1 = 128 turns, tap 2 = 128+23 = 151 turns.
+    name : str
+        Optional label (e.g. "primary", "secondary").
+    """
+    awg: "int | list[int]"
+    taps: list[int]
+    name: str = "winding"
+
+    def __post_init__(self):
+        if len(self.taps) < 2:
+            raise ValueError("taps must have at least 2 elements (first is always 0)")
+        if self.taps[0] != 0:
+            raise ValueError("First element of taps must be 0")
+        if any(t < 0 for t in self.taps):
+            raise ValueError("All tap increments must be >= 0")
+        n_segments = len(self.taps) - 1
+        if isinstance(self.awg, int):
+            self.awg = [self.awg] * n_segments
+        if len(self.awg) != n_segments:
+            raise ValueError(
+                f"awg list length ({len(self.awg)}) must match number of segments ({n_segments})"
+            )
+
+    @property
+    def total_turns(self) -> int:
+        """Total turns for the highest tap."""
+        return sum(self.taps)
+
+    def cumulative_turns(self) -> list[int]:
+        """Cumulative turns at each tap point (1-indexed tap numbers)."""
+        result = []
+        total = 0
+        for t in self.taps[1:]:
+            total += t
+            result.append(total)
+        return result
+
+    def segment_turns(self) -> list[int]:
+        """Turns in each physical segment (same as taps[1:])."""
+        return list(self.taps[1:])
+
+
+@dataclass
+class LayerResult:
+    """Results for one winding layer."""
+    layer_index: int       # 0-based
+    turns_capacity: int    # max turns that fit in this layer
+    turns_used: int        # turns actually placed in this layer
+    L_turn_mm: float       # mean turn length for this layer (mm)
+    wire_length_mm: float  # total wire length for turns in this layer (mm)
+
+
+@dataclass
+class SegmentResult:
+    """Results for one tap segment."""
+    segment_index: int     # 0-based
+    awg: int
+    turns: int
+    layers: list[LayerResult]
+    wire_length_m: float
+    resistance_ohm: float
+    weight_g: float
+    fits: bool             # True if all turns fit within fill factor and geometry
+
+
+@dataclass
+class WindingResult:
+    """Full winding analysis result."""
+    spec: WindingSpec
+    core: ToroidCore
+    fill_factor: float
+    max_turns_fill: int          # turns allowed by fill factor
+    max_turns_geometry: int      # turns allowed by geometry alone
+    segments: list[SegmentResult]
+    total_wire_length_m: float
+    total_resistance_ohm: float
+    total_weight_g: float
+    feasible: bool               # True if all segments fit
+    # Cost fields (None when AWG not in cost table)
+    spools_required: Optional[int] = None
+    cost_usd: Optional[float] = None
+    # Internal winding position at end of this winding (for chaining)
+    _end_layer: int = 0
+    _end_turns_in_layer: int = 0
+    _end_total_turns: int = 0
+
+    def summary(self) -> str:
+        # Header: show AWG as single value if uniform, else show range
+        awg_list = self.spec.awg  # always a list after __post_init__
+        awg_str = str(awg_list[0]) if len(set(awg_list)) == 1 else str(awg_list)
+        lines = [
+            f"=== Winding: {self.spec.name} (AWG {awg_str}) ===",
+            f"  Core: ID={self.core.ID_mm}mm  OD={self.core.OD_mm}mm  H={self.core.height_mm}mm",
+            f"  Fill factor: {self.fill_factor*100:.0f}%  "
+            f"(max turns by fill={self.max_turns_fill}, by geometry={self.max_turns_geometry})",
+            f"  Total turns: {self.spec.total_turns}  |  Feasible: {'YES' if self.feasible else 'NO -- DOES NOT FIT'}",
+            "",
+        ]
+        cumulative = 0
+        for seg in self.segments:
+            cumulative += seg.turns
+            tap_num = seg.segment_index + 1
+            status = "OK" if seg.fits else "OVERFLOW"
+            # Show AWG per tap only when mixed
+            awg_tag = f"  AWG {seg.awg}" if len(set(awg_list)) > 1 else ""
+            lines.append(
+                f"  Tap {tap_num}: {seg.turns:>4} turns  "
+                f"(cumulative {cumulative:>4}){awg_tag}  "
+                f"{seg.wire_length_m:7.2f} m  "
+                f"{seg.resistance_ohm*1000:8.3f} mOhm  "
+                f"{seg.weight_g:7.3f} g  [{status}]"
+            )
+            for lr in seg.layers:
+                lines.append(
+                    f"    layer {lr.layer_index:2d}: capacity={lr.turns_capacity:3d}  "
+                    f"used={lr.turns_used:3d}  "
+                    f"MTL={lr.L_turn_mm:.1f}mm  "
+                    f"wire={lr.wire_length_mm/1000:.3f}m"
+                )
+        cost_str = ""
+        if self.cost_usd is not None:
+            cost_str = f"  |  ${self.cost_usd:.2f}"
+        lines += [
+            "",
+            f"  TOTAL: {self.total_wire_length_m:.2f} m  |  "
+            f"{self.total_resistance_ohm*1000:.3f} mOhm  |  "
+            f"{self.total_weight_g:.2f} g"
+            + cost_str,
+        ]
+        return "\n".join(lines)
+
+
+def _layer_geometry(core: ToroidCore, wire: WireSpec, layer_index: int
+                    ) -> tuple[int, float]:
+    """
+    Calculate the capacity and mean turn length for a given layer index.
+
+    Returns
+    -------
+    (turns_capacity, L_turn_mm)
+        turns_capacity : maximum turns that fit in this layer (0 if no room)
+        L_turn_mm      : mean path length of one turn in this layer (mm)
+    """
+    n = layer_index
+    d = wire.diameter_mm
+
+    # Wire centres sit at ID/2 - (n+0.5)*d from the bore axis
+    wire_centre_radius = core.ID_mm / 2.0 - (n + 0.5) * d
+    # Available height for this layer (each layer adds one wire diameter
+    # of build-up on top and bottom faces, reducing the straight section)
+    avail_height = core.height_mm - 2 * n * d
+
+    if wire_centre_radius <= 0 or avail_height <= 0:
+        return 0, 0.0
+
+    # Turns per layer = how many wire diameters fit around the inner
+    # circumference at the wire centre radius.  The height only constrains
+    # how many *layers* can be wound, not how many turns per layer.
+    inner_circumference = 2 * math.pi * wire_centre_radius
+    turns_capacity = int(inner_circumference / d)   # floor
+
+    if turns_capacity <= 0:
+        return 0, 0.0
+
+    # Mean turn length: inner + outer semicircle + two straight sections
+    outer_radius = core.OD_mm / 2.0 + (n + 0.5) * d
+    L_turn_mm = math.pi * (wire_centre_radius + outer_radius) + 2 * avail_height
+
+    return turns_capacity, L_turn_mm
+
+
+def analyse_winding(
+    core: ToroidCore,
+    spec: WindingSpec,
+    fill_factor: float = 0.35,
+    start_layer: int = 0,
+    start_turns_in_layer: int = 0,
+    start_total_turns: int = 0,
+) -> WindingResult:
+    """
+    Analyse feasibility and compute wire parameters for a winding.
+
+    Parameters
+    ----------
+    core : ToroidCore
+    spec : WindingSpec
+    fill_factor : float
+        Maximum fraction of the bore window area that may be filled.
+        Default 0.35 (35%).
+    start_layer : int
+        Layer index to begin winding on.  Pass the previous winding's
+        _end_layer to continue mid-layer across windings.
+    start_turns_in_layer : int
+        Turns already consumed in start_layer from a previous winding.
+    start_total_turns : int
+        Cumulative turns already placed (used for fill-factor accounting
+        across windings).
+
+    Returns
+    -------
+    WindingResult
+    """
+    # Resolve per-segment wire specs (spec.awg is always a list after __post_init__)
+    seg_wires = [WireSpec.from_awg(a) for a in spec.awg]
+
+    # For informational max_turns_fill / max_turns_geometry use the thickest
+    # wire (largest diameter) as a conservative lower bound.
+    wire_ref = max(seg_wires, key=lambda w: w.diameter_mm)
+    wire_area_ref = math.pi * (wire_ref.diameter_mm / 2.0) ** 2
+    max_turns_fill = int(fill_factor * core.window_area_mm2 / wire_area_ref)
+
+    max_turns_geometry = 0
+    layer_idx = 0
+    while True:
+        cap, _ = _layer_geometry(core, wire_ref, layer_idx)
+        if cap == 0:
+            break
+        max_turns_geometry += cap
+        layer_idx += 1
+
+    effective_max = max_turns_fill  # fill factor is the binding constraint
+
+    segments: list[SegmentResult] = []
+    overall_feasible = True
+
+    # Track position across segments: which layer we're on and how many
+    # turns have been placed in the current layer already.
+    current_layer = start_layer
+    turns_in_current_layer = start_turns_in_layer
+    total_turns_placed = start_total_turns
+    current_layer_wire_diameter: Optional[float] = None   # set on first use
+
+    for seg_idx, (seg_turns, wire) in enumerate(zip(spec.segment_turns(), seg_wires)):
+        seg_layers: list[LayerResult] = []
+        seg_wire_length_mm = 0.0
+        seg_fits = True
+        turns_remaining = seg_turns
+
+        while turns_remaining > 0:
+            # If the wire gauge changed from what's already on this layer,
+            # start a fresh layer so gauges are never mixed on one layer.
+            if (current_layer_wire_diameter is not None
+                    and wire.diameter_mm != current_layer_wire_diameter
+                    and turns_in_current_layer > 0):
+                current_layer += 1
+                turns_in_current_layer = 0
+
+            cap, L_turn = _layer_geometry(core, wire, current_layer)
+            current_layer_wire_diameter = wire.diameter_mm
+
+            if cap == 0:
+                # No more geometry room at all
+                seg_fits = False
+                overall_feasible = False
+                # Record overflow as a single "layer" with 0 capacity
+                seg_layers.append(LayerResult(
+                    layer_index=current_layer,
+                    turns_capacity=0,
+                    turns_used=turns_remaining,
+                    L_turn_mm=0.0,
+                    wire_length_mm=0.0,
+                ))
+                turns_remaining = 0
+                break
+
+            available_in_layer = cap - turns_in_current_layer
+            if available_in_layer <= 0:
+                # Layer is full, advance
+                current_layer += 1
+                turns_in_current_layer = 0
+                continue
+
+            # Check fill-factor cap
+            if total_turns_placed >= effective_max:
+                seg_fits = False
+                overall_feasible = False
+                seg_layers.append(LayerResult(
+                    layer_index=current_layer,
+                    turns_capacity=available_in_layer,
+                    turns_used=turns_remaining,
+                    L_turn_mm=L_turn,
+                    wire_length_mm=0.0,
+                ))
+                turns_remaining = 0
+                break
+
+            turns_to_place = min(turns_remaining,
+                                 available_in_layer,
+                                 effective_max - total_turns_placed)
+            wire_len = turns_to_place * L_turn
+
+            seg_layers.append(LayerResult(
+                layer_index=current_layer,
+                turns_capacity=cap,
+                turns_used=turns_to_place,
+                L_turn_mm=L_turn,
+                wire_length_mm=wire_len,
+            ))
+
+            seg_wire_length_mm += wire_len
+            turns_remaining -= turns_to_place
+            turns_in_current_layer += turns_to_place
+            total_turns_placed += turns_to_place
+
+            if turns_in_current_layer >= cap:
+                current_layer += 1
+                turns_in_current_layer = 0
+
+        seg_wire_length_m = seg_wire_length_mm / 1000.0
+        seg_resistance = (wire.ohm_per_km / 1000.0) * seg_wire_length_m
+        # Weight from conductor area × length × density
+        seg_weight_g = wire.area_mm2 * seg_wire_length_mm * _COPPER_DENSITY_G_MM3
+
+        segments.append(SegmentResult(
+            segment_index=seg_idx,
+            awg=wire.awg,
+            turns=seg_turns,
+            layers=seg_layers,
+            wire_length_m=seg_wire_length_m,
+            resistance_ohm=seg_resistance,
+            weight_g=seg_weight_g,
+            fits=seg_fits,
+        ))
+
+    total_wire_m = sum(s.wire_length_m for s in segments)
+    total_resistance = sum(s.resistance_ohm for s in segments)
+    total_weight = sum(s.weight_g for s in segments)
+
+    # Cost: pro-rated per segment (each may have different AWG/price)
+    spools_required = None
+    cost_usd = None
+    seg_costs = []
+    for seg in segments:
+        w = WireSpec.from_awg(seg.awg)
+        if w.cost_per_m is not None:
+            seg_costs.append(w.cost_per_m * seg.wire_length_m)
+    if seg_costs:
+        cost_usd = sum(seg_costs)
+
+    return WindingResult(
+        spec=spec,
+        core=core,
+        fill_factor=fill_factor,
+        max_turns_fill=max_turns_fill,
+        max_turns_geometry=max_turns_geometry,
+        segments=segments,
+        total_wire_length_m=total_wire_m,
+        total_resistance_ohm=total_resistance,
+        total_weight_g=total_weight,
+        feasible=overall_feasible,
+        spools_required=spools_required,
+        cost_usd=cost_usd,
+        _end_layer=current_layer,
+        _end_turns_in_layer=turns_in_current_layer,
+        _end_total_turns=total_turns_placed,
+    )
+
+
+def design_transformer(
+    core: ToroidCore,
+    windings: list[WindingSpec],
+    fill_factor: float = 0.35,
+) -> list[WindingResult]:
+    """
+    Analyse multiple windings on the same core, packing them continuously
+    mid-layer.  Each winding picks up exactly where the previous one left off.
+
+    The fill_factor budget is shared across all windings: once the total turns
+    placed (regardless of AWG) reaches fill_factor * window_area / wire_area,
+    no more turns can be added.  Because different AWGs have different wire
+    areas, the budget is tracked in terms of area consumed.
+
+    Parameters
+    ----------
+    core : ToroidCore
+    windings : list[WindingSpec]
+        Windings in winding order.
+    fill_factor : float
+        Maximum fraction of bore window area consumed by ALL windings combined.
+
+    Returns
+    -------
+    list[WindingResult]  (one per winding, in order)
+    """
+    results: list[WindingResult] = []
+    cur_layer = 0
+    cur_turns_in_layer = 0
+    # Track consumed window area (mm²) across all windings to enforce the
+    # shared fill-factor budget regardless of AWG.
+    consumed_area_mm2 = 0.0
+    budget_area_mm2 = fill_factor * core.window_area_mm2
+
+    for spec in windings:
+        # Use the thickest wire in this spec for area-budget accounting
+        seg_wires = [WireSpec.from_awg(a) for a in spec.awg]
+        wire_ref = max(seg_wires, key=lambda w: w.diameter_mm)
+        wire_area_mm2 = math.pi * (wire_ref.diameter_mm / 2.0) ** 2
+        # Express the remaining area budget as an equivalent turn count for
+        # this winding's reference wire gauge, then back-calculate the
+        # synthetic start offset so that analyse_winding's fill-factor
+        # headroom check is correct.
+        turns_already_equivalent = int(consumed_area_mm2 / wire_area_mm2)
+
+        result = analyse_winding(
+            core=core,
+            spec=spec,
+            fill_factor=fill_factor,
+            start_layer=cur_layer,
+            start_turns_in_layer=cur_turns_in_layer,
+            start_total_turns=turns_already_equivalent,
+        )
+        results.append(result)
+
+        # Advance shared layer position
+        cur_layer = result._end_layer
+        cur_turns_in_layer = result._end_turns_in_layer
+        # Update consumed area using actual turns placed in this winding
+        turns_placed = result._end_total_turns - turns_already_equivalent
+        consumed_area_mm2 = min(budget_area_mm2,
+                                consumed_area_mm2 + turns_placed * wire_area_mm2)
+
+    return results
+
+
+# ---------------------------------------------------------------------------
+# Example / smoke-test
+# ---------------------------------------------------------------------------
+
+if __name__ == "__main__":
+    core = ToroidCore(ID_mm=21.5, OD_mm=46.5, height_mm=22.8)
+
+    primary = WindingSpec(
+        awg=[22, 22],      # tap 1 = AWG 22, tap 2 = AWG 24
+        taps=[0, 25, 50],
+        name="primary",
+    )
+
+    secondary = WindingSpec(
+        awg=[22, 22, 22, 26],            # uniform gauge
+        taps=[0, 100, 50, 50, 50],
+        name="secondary",
+    )
+
+    results = design_transformer(core, [primary, secondary], fill_factor=0.35)
+    for result in results:
+        print(result.summary())
+        print()

draw_toroid.py

@@ -0,0 +1,368 @@
+"""
+draw_toroid.py
+==============
+Top-down (axial) cross-section drawing of a toroidal transformer.
+
+The view is looking straight down the axis of the toroid.
+
+What you see in this view
+-------------------------
+  - Grey annulus  : the core material (ID to OD)
+  - White disc    : the bore window (diameter = ID)
+  - Coloured dots : wire cross-sections packed against the inner bore wall,
+                    one concentric ring per winding layer
+
+Physical basis: each wire turn passes through the bore.  Looking down the
+axis you see the wire cross-sections arranged in rings inside the bore.
+Layer 0 (first layer) sits closest to the bore wall; successive layers sit
+inward.  Wire centres for layer n are at:
+
+    r(n) = ID/2 - (n + 0.5) * d_wire
+
+Turns are spaced by d_wire around the circumference of that ring.  Up to
+turns_capacity dots are drawn per ring (the physical fill of that layer).
+
+Usage
+-----
+    from designer import ToroidCore, WindingSpec, design_transformer
+    from draw_toroid import draw_toroid
+
+    core    = ToroidCore(ID_mm=21.5, OD_mm=46.5, height_mm=22.8)
+    primary = WindingSpec(awg=22, taps=[0, 25, 50], name="primary")
+    sec     = WindingSpec(awg=22, taps=[0, 100, 50], name="secondary")
+    results = design_transformer(core, [primary, secondary])
+
+    draw_toroid(core, results, "toroid.png")
+"""
+
+from __future__ import annotations
+
+import math
+from typing import List
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import matplotlib.patches as mpatches
+from matplotlib.patches import Circle
+import numpy as np
+
+from designer import ToroidCore, WindingResult, WireSpec
+
+
+# ---------------------------------------------------------------------------
+# Colour scheme  (one per winding)
+# ---------------------------------------------------------------------------
+_WINDING_COLORS = [
+    "#1565C0",   # dark blue  — winding 0
+    "#B71C1C",   # dark red   — winding 1
+    "#2E7D32",   # dark green — winding 2
+    "#E65100",   # dark orange— winding 3
+    "#6A1B9A",   # dark purple— winding 4
+    "#00838F",   # dark cyan  — winding 5
+]
+
+_CORE_FACE_COLOR = "#9E9E9E"   # core material
+_CORE_EDGE_COLOR = "#212121"
+_WIRE_EDGE_COLOR = "#000000"
+_BG_COLOR        = "white"
+
+
+# ---------------------------------------------------------------------------
+# Public API
+# ---------------------------------------------------------------------------
+
+def draw_toroid(
+    core: ToroidCore,
+    winding_results: List[WindingResult],
+    output_path: str = "toroid.png",
+    dpi: int = 180,
+    fig_size_mm: float = 180.0,
+) -> None:
+    """
+    Render a to-scale top-down cross-section of the transformer as a PNG.
+
+    Parameters
+    ----------
+    core : ToroidCore
+    winding_results : list[WindingResult]
+        Output of design_transformer() or analyse_winding().
+    output_path : str
+        Destination file (PNG).
+    dpi : int
+        Render resolution.
+    fig_size_mm : float
+        Square canvas size in mm.
+    """
+    OD = core.OD_mm
+    ID = core.ID_mm
+
+    # Canvas extent: just enough margin around the OD
+    margin = OD * 0.12
+    extent = OD / 2.0 + margin      # mm from centre to canvas edge
+
+    fig_in = fig_size_mm / 25.4
+    fig, ax = plt.subplots(figsize=(fig_in, fig_in), facecolor=_BG_COLOR)
+    ax.set_facecolor(_BG_COLOR)
+    ax.set_aspect("equal")
+    ax.set_xlim(-extent, extent)
+    ax.set_ylim(-extent, extent)
+    ax.set_xlabel("mm", fontsize=8)
+    ax.set_ylabel("mm", fontsize=8)
+    ax.tick_params(labelsize=7)
+    ax.set_title(
+        f"Toroid cross-section  "
+        f"ID={ID:.1f} mm  OD={OD:.1f} mm  H={core.height_mm:.1f} mm",
+        fontsize=9, pad=6,
+    )
+
+    # -----------------------------------------------------------------------
+    # Core annulus
+    # -----------------------------------------------------------------------
+    core_ring = mpatches.Annulus(
+        xy=(0, 0),
+        r=OD / 2.0,
+        width=(OD - ID) / 2.0,
+        facecolor=_CORE_FACE_COLOR,
+        edgecolor=_CORE_EDGE_COLOR,
+        linewidth=1.2,
+        zorder=2,
+        label=f"Core  ID={ID:.1f} OD={OD:.1f} H={core.height_mm:.1f} mm",
+    )
+    ax.add_patch(core_ring)
+
+    # Bore window
+    bore = Circle(
+        (0, 0), ID / 2.0,
+        facecolor=_BG_COLOR,
+        edgecolor=_CORE_EDGE_COLOR,
+        linewidth=0.8,
+        zorder=3,
+    )
+    ax.add_patch(bore)
+
+
+    # -----------------------------------------------------------------------
+    # Build layer_index -> wire_centre_radius map (accounts for mixed gauges)
+    #
+    # Layers are wound in order from the bore wall inward.  Each layer
+    # consumes its own wire diameter of radial space.  We walk through all
+    # layers across all windings in ascending layer_index order and accumulate
+    # the actual radial depth so the centre radius is correct regardless of
+    # whether wire gauge changes between segments or windings.
+    # -----------------------------------------------------------------------
+
+    # Collect (layer_index, wire_diameter) for every active layer, all windings
+    _layer_d: dict[int, float] = {}
+    for wr in winding_results:
+        for seg in wr.segments:
+            wire = WireSpec.from_awg(seg.awg)
+            for lr in seg.layers:
+                if lr.turns_used > 0 and lr.turns_capacity > 0:
+                    # First winding to touch a layer defines its wire gauge
+                    if lr.layer_index not in _layer_d:
+                        _layer_d[lr.layer_index] = wire.diameter_mm
+
+    # Walk layers in order to accumulate radial depth from the bore wall
+    layer_centre_r: dict[int, float] = {}
+    depth = 0.0
+    for n in sorted(_layer_d):
+        d = _layer_d[n]
+        layer_centre_r[n] = ID / 2.0 - depth - d / 2.0
+        depth += d
+
+    # -----------------------------------------------------------------------
+    # Pre-compute per-layer start angles (shared across all windings) so
+    # the total winding arc on each layer is centred at the top (π/2).
+    # -----------------------------------------------------------------------
+    _layer_total_used: dict[int, int] = {}
+    for _wr in winding_results:
+        for _seg in _wr.segments:
+            for _lr in _seg.layers:
+                if _lr.turns_used > 0 and _lr.turns_capacity > 0:
+                    n_ = _lr.layer_index
+                    _layer_total_used[n_] = _layer_total_used.get(n_, 0) + _lr.turns_used
+
+    # Step = d/r (touching circles), start at 0, wind continuously.
+    # All segments share the same running angle per layer — they just
+    # pick up where the previous segment left off and keep going around.
+    layer_angle_step: dict[int, float] = {}
+    layer_next_angle: dict[int, float] = {}
+    for n_ in _layer_total_used:
+        r_ = layer_centre_r.get(n_, 0.0)
+        d_ = _layer_d.get(n_, 1.0)
+        if r_ > 0:
+            layer_angle_step[n_] = d_ / r_
+            layer_next_angle[n_] = 0.0
+
+    legend_handles: list[mpatches.Patch] = []
+
+    for w_idx, wr in enumerate(winding_results):
+        base_color = _WINDING_COLORS[w_idx % len(_WINDING_COLORS)]
+        awg_list = wr.spec.awg
+        awg_str = str(awg_list[0]) if len(set(awg_list)) == 1 else f"{awg_list[0]}-{awg_list[-1]}"
+        n_segs = len(wr.spec.taps) - 1
+        feasible_str = "OK" if wr.feasible else "DOES NOT FIT"
+        legend_handles.append(mpatches.Patch(
+            facecolor=base_color,
+            edgecolor=_WIRE_EDGE_COLOR,
+            linewidth=0.8,
+            label=f"{wr.spec.name}  AWG {awg_str}  {n_segs} seg(s)  "
+                  f"{wr.spec.total_turns} turns  [{feasible_str}]",
+        ))
+
+        # Colour shading: outermost layers → light, innermost → dark.
+        for seg in wr.segments:
+            wire = WireSpec.from_awg(seg.awg)
+            d = wire.diameter_mm
+
+            # Shade by segment index: seg 0 → lightest, last seg → base colour.
+            # Use full 0..1 range across all segments for maximum separation.
+            frac = seg.segment_index / max(1, n_segs - 1)
+            seg_color = _lighten(base_color, 0.55 * (1.0 - frac))
+
+            for lr in seg.layers:
+                if lr.turns_used == 0 or lr.turns_capacity == 0:
+                    continue
+
+                n = lr.layer_index
+                r = layer_centre_r.get(n, 0.0)
+                if r <= 0:
+                    continue
+
+                angle_step = layer_angle_step.get(n, 2.0 * math.pi / max(lr.turns_used, 1))
+                start_angle = layer_next_angle.get(n, 0.0)
+
+                # Draw at true wire radius. Circles are evenly spaced over
+                # 360° so they tile the ring; they may overlap on outer
+                # layers (where arc spacing > wire diameter) or underlap on
+                # very dense inner layers, but the count and colour are correct.
+                draw_r = d / 2.0
+
+                for i in range(lr.turns_used):
+                    a = start_angle + i * angle_step
+                    x = r * math.cos(a)
+                    y = r * math.sin(a)
+                    ax.add_patch(Circle(
+                        (x, y), draw_r,
+                        facecolor=seg_color,
+                        edgecolor=_WIRE_EDGE_COLOR,
+                        linewidth=0.35,
+                        alpha=0.90,
+                        zorder=10 + n,
+                    ))
+
+                # Advance the angle for the next segment on this layer
+                layer_next_angle[n] = start_angle + lr.turns_used * angle_step
+
+            # Segment legend entry
+            awg_tag = f"  AWG {seg.awg}" if len(set(awg_list)) > 1 else ""
+            n_active_layers = sum(
+                1 for lr in seg.layers if lr.turns_used > 0 and lr.turns_capacity > 0
+            )
+            legend_handles.append(mpatches.Patch(
+                facecolor=seg_color,
+                edgecolor=_WIRE_EDGE_COLOR,
+                linewidth=0.6,
+                label=(f"  tap {seg.segment_index + 1}: {seg.turns} turns{awg_tag}"
+                       f"  {n_active_layers} layer(s)"
+                       f"  {seg.resistance_ohm * 1e3:.1f} mOhm"),
+            ))
+
+    # -----------------------------------------------------------------------
+    # Scale bar (bottom-centre)
+    # -----------------------------------------------------------------------
+    bar_mm = _nice_bar(OD)
+    bx = -bar_mm / 2
+    by = -(extent * 0.87)
+    tick = extent * 0.012
+    ax.plot([bx, bx + bar_mm], [by, by], color="black", lw=1.5, zorder=20)
+    ax.plot([bx,          bx],          [by - tick, by + tick], color="black", lw=1.5, zorder=20)
+    ax.plot([bx + bar_mm, bx + bar_mm], [by - tick, by + tick], color="black", lw=1.5, zorder=20)
+    ax.text(bx + bar_mm / 2, by - tick * 1.5, f"{bar_mm:.0f} mm",
+            ha="center", va="top", fontsize=7)
+
+    # -----------------------------------------------------------------------
+    # Legend (right of axes)
+    # -----------------------------------------------------------------------
+    ax.legend(
+        handles=legend_handles,
+        loc="upper left",
+        bbox_to_anchor=(1.02, 1.0),
+        borderaxespad=0,
+        fontsize=7,
+        frameon=True,
+        title="Windings & taps",
+        title_fontsize=8,
+    )
+
+    fig.savefig(output_path, dpi=dpi, bbox_inches="tight", facecolor=_BG_COLOR)
+    plt.close(fig)
+    print(f"Saved: {output_path}")
+
+
+# ---------------------------------------------------------------------------
+# Geometry helpers
+# ---------------------------------------------------------------------------
+
+def _ring_positions(radius: float, n: int,
+                    d_wire: float = 0.0) -> list[tuple[float, float]]:
+    """
+    Place n wire centres uniformly around a full 360° circle.
+
+    n should equal turns_capacity (floor(2πr/d)), so the spacing between
+    adjacent centres is approximately d_wire — the circles are tangent with
+    no visible gap.  Always covers the full circumference.
+    """
+    if n <= 0 or radius <= 0:
+        return []
+    angles = np.linspace(0.0, 2 * math.pi, n, endpoint=False)
+    return [(radius * math.cos(a), radius * math.sin(a)) for a in angles]
+
+
+
+def _nice_bar(OD_mm: float) -> float:
+    """Round scale-bar length ≈ 20% of OD."""
+    raw = OD_mm * 0.2
+    for v in [1, 2, 5, 10, 20, 25, 50]:
+        if v >= raw:
+            return float(v)
+    return float(round(raw / 10) * 10)
+
+
+def _lighten(hex_color: str, amount: float) -> str:
+    """
+    Mix hex_color towards white by `amount` (0=unchanged, 1=white).
+    Used to shade tap segments within a winding.
+    """
+    amount = max(0.0, min(1.0, amount))
+    h = hex_color.lstrip("#")
+    r, g, b = (int(h[i:i+2], 16) / 255.0 for i in (0, 2, 4))
+    r = r + (1 - r) * amount
+    g = g + (1 - g) * amount
+    b = b + (1 - b) * amount
+    return "#{:02X}{:02X}{:02X}".format(int(r * 255), int(g * 255), int(b * 255))
+
+
+# ---------------------------------------------------------------------------
+# Smoke-test
+# ---------------------------------------------------------------------------
+
+if __name__ == "__main__":
+    from designer import WindingSpec, design_transformer
+
+    core = ToroidCore(ID_mm=21.5, OD_mm=46.5, height_mm=22.8)
+
+    primary = WindingSpec(
+        awg=[22, 22],
+        taps=[0, 25, 50],
+        name="primary",
+    )
+    secondary = WindingSpec(
+        awg=[22, 22, 22, 26],
+        taps=[0, 100, 50, 50, 50],
+        name="secondary",
+    )
+
+    results = design_transformer(core, [primary, secondary])
+    draw_toroid(core, results, "toroid.png", dpi=180)

model.py

@@ -1,5 +1,6 @@
 from dataclasses import dataclass, field
 from typing import List, Tuple, Optional, Dict
+from core import core_loss_Pv
 
 
 @dataclass
@@ -18,6 +19,9 @@ class TransformerModel:
     # Turn counts (total)
     Np_total: int
     Ns_total: int
+    # Core loss model parameters (optional)
+    Ve_mm3: float = 0.0          # effective core volume in mm^3 (for core loss calculation)
+    use_core_loss_model: bool = False  # if True, calculate core loss from B and f
     # 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])
@@ -173,8 +177,26 @@ class TransformerModel:
         # Output power (resistive load)
         P_out = Vs_rms_ideal**2 / load_ohms
 
+        # Core loss calculation
+        if self.use_core_loss_model and self.Ve_mm3 > 0:
+            # Calculate core loss from B and f using the model
+            # Convert volume to m^3
+            Ve_m3 = self.Ve_mm3 * 1e-9  # mm^3 to m^3
+
+            # Get core loss density in kW/m^3
+            try:
+                Pv_kW_m3 = core_loss_Pv(B_peak_T, freq_hz)
+                # Convert to watts
+                core_loss_calculated = Pv_kW_m3 * Ve_m3 * 1000  # kW/m^3 * m^3 * 1000 = W
+            except (ValueError, Exception) as e:
+                # If core loss model fails (e.g., B or f out of range), fall back to provided value
+                core_loss_calculated = core_loss_W
+        else:
+            # Use provided core loss value
+            core_loss_calculated = core_loss_W
+
         # Total input power (approx)
-        P_in = P_out + P_cu_total + core_loss_W
+        P_in = P_out + P_cu_total + core_loss_calculated
 
         efficiency = P_out / P_in if P_in > 0 else 1.0
 
@@ -191,7 +213,7 @@ class TransformerModel:
             "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_core_W": core_loss_calculated,
             "P_in_W": P_in,
             "efficiency": efficiency,
             "primary_tap": primary_tap,

optimizer.py

@@ -57,6 +57,7 @@ class TransformerOptimizer:
         Vp_step: float = 0.5,
         B_max_T: float = 0.3,
         Vs_max: float = float('inf'),
+        Is_max: float = float('inf'),
         core_loss_W: float = 0.0,
         power_tolerance_percent: float = 2.0,
         fallback_max_power: bool = True,
@@ -71,6 +72,7 @@ class TransformerOptimizer:
         - 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)
+        - Is_max: Maximum allowed secondary current (A)
         - 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
@@ -115,6 +117,7 @@ class TransformerOptimizer:
                 Vp_step=Vp_step,
                 B_max_T=B_max_T,
                 Vs_max=Vs_max,
+                Is_max=Is_max,
                 core_loss_W=core_loss_W,
                 power_tolerance_percent=power_tolerance_percent,
                 fallback_max_power=fallback_max_power,
@@ -151,6 +154,7 @@ class TransformerOptimizer:
         Vp_step: float,
         B_max_T: float,
         Vs_max: float,
+        Is_max: float,
         core_loss_W: float,
         power_tolerance_percent: float,
         fallback_max_power: bool,
@@ -195,6 +199,9 @@ class TransformerOptimizer:
             if sim["Vs_rms"] > Vs_max:
                 continue
 
+            if sim["Is_rms"] > Is_max:
+                continue
+
             power_error = abs(sim["P_out_W"] - target_power_W) / target_power_W * 100
 
             result = OptimizationResult(

requirements.txt

@@ -0,0 +1,2 @@
+flask>=2.3.0
+numpy>=1.24.0

sim.py

@@ -86,9 +86,9 @@ 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,
+    load_ohms=10.0,
-    target_power_W=10.0,
+    target_power_W=25.0,
-    freq_hz=45000.0,
+    freq_hz=256.0,
     Vp_min=1.0,
     Vp_max=36.0,
     Vp_step=0.5,

sim_toroid.py

@@ -0,0 +1,970 @@
+"""
+sim_toroid.py
+=============
+Simulate a toroidal transformer designed with designer.py.
+
+Given a pair of WindingResult objects (primary / secondary) produced by
+designer.py, this module lets you evaluate any operating point and check
+it against a set of constraints.
+
+Typical usage
+-------------
+    from designer import ToroidCore, WindingSpec, design_transformer
+    from sim_toroid import ToroidSimulator, SimConstraints
+
+    core = ToroidCore(
+        ID_mm=21.5, OD_mm=46.5, height_mm=22.8,
+        Ae_mm2=297.0,           # effective cross-section (mm²)
+        Ve_mm3=18500.0,         # effective volume (mm³)
+        Pv_func=my_pv_func,     # Pv_func(f_hz, B_T) -> kW/m³  (or None)
+    )
+    primary  = WindingSpec(awg=22, taps=[0, 25, 50], name="primary")
+    secondary = WindingSpec(awg=22, taps=[0, 100, 50], name="secondary")
+    pri_result, sec_result = design_transformer(core, [primary, secondary])
+
+    sim = ToroidSimulator(core=core, primary=pri_result, secondary=sec_result)
+
+    constraints = SimConstraints(
+        B_max_T=0.3,
+        Vp_max=50.0,
+        Vs_max=30.0,
+        Ip_max=5.0,
+        Is_max=3.0,
+        P_out_max_W=50.0,
+    )
+
+    result = sim.simulate(
+        Vp_rms=24.0,
+        freq_hz=50_000.0,
+        primary_tap=2,
+        secondary_tap=1,
+        Z_load=(10.0, 0.0),          # 10 Ω purely resistive
+        constraints=constraints,
+    )
+    print(result)
+
+Physics model
+-------------
+The model uses the standard ideal-transformer approximation with
+first-order winding resistance correction:
+
+  * Flux density:   B_peak = Vp / (4.44 · Np · Ae · f)
+  * Ideal secondary voltage: Vs_ideal = Vp · (Ns / Np)
+  * Secondary circuit:  Vs_oc drives (Rs + Z_load) series circuit
+        Is  = Vs_ideal / (Z_load + Rs)
+        Vs  = Is · Z_load
+  * Primary current reflected from secondary plus magnetising (ignored):
+        Ip  = Is · (Ns / Np)          [for turns-ratio correction]
+        Then added primary drop:
+        Vp_actual ≈ Vp  (Vp is the applied voltage, copper drop on primary
+                         reduces flux slightly; accounted by reducing Vp seen
+                         by the core by Ip·Rp)
+  * Core loss computed via core.Pv_func(f, B) [kW/m³] · core.Ve_mm3 [mm³]
+"""
+
+from __future__ import annotations
+
+import cmath
+import math
+from dataclasses import dataclass, field
+from typing import Callable, ClassVar, Optional, Tuple
+
+from designer import ToroidCore, WindingResult
+
+
+# ---------------------------------------------------------------------------
+# Result dataclass
+# ---------------------------------------------------------------------------
+
+@dataclass
+class SimResult:
+    """Operating-point results from ToroidSimulator.simulate()."""
+
+    # Inputs (echoed)
+    Vp_rms: float
+    freq_hz: float
+    primary_tap: int
+    secondary_tap: int
+
+    # Turns
+    Np_eff: int
+    Ns_eff: int
+    turns_ratio: float          # Ns / Np
+
+    # Magnetics
+    B_peak_T: float             # peak flux density in core
+
+    # Voltages (RMS magnitudes)
+    Vs_rms: float               # voltage across load
+    Vp_rms_applied: float       # = Vp_rms (the applied primary voltage)
+
+    # Currents (RMS magnitudes)
+    Ip_rms: float               # primary current
+    Is_rms: float               # secondary current
+
+    # Phase angle of load current relative to secondary voltage (radians)
+    load_phase_rad: float
+
+    # Winding resistances (ohms)
+    Rp: float                   # primary winding resistance for active turns
+    Rs: float                   # secondary winding resistance for active turns
+
+    # Power breakdown (watts)
+    P_out_W: float              # real power delivered to load
+    P_cu_primary_W: float       # primary copper loss
+    P_cu_secondary_W: float     # secondary copper loss
+    P_cu_W: float               # total copper loss
+    P_core_W: float             # core loss
+    P_in_W: float               # total input power
+    efficiency: float           # P_out / P_in  (0–1)
+
+    # Constraint violations (empty list = all OK)
+    violations: list[str] = field(default_factory=list)
+
+    @property
+    def feasible(self) -> bool:
+        return len(self.violations) == 0
+
+    def __str__(self) -> str:
+        lines = [
+            f"=== Simulation result (tap P{self.primary_tap}/S{self.secondary_tap}, "
+            f"{self.freq_hz/1e3:.1f} kHz, Vp={self.Vp_rms:.2f} V) ===",
+            f"  Turns          : Np={self.Np_eff}  Ns={self.Ns_eff}  ratio={self.turns_ratio:.4f}",
+            f"  Flux density   : B_peak = {self.B_peak_T*1e3:.2f} mT",
+            f"  Winding R      : Rp = {self.Rp*1e3:.3f} mOhm   Rs = {self.Rs*1e3:.3f} mOhm",
+            f"  Voltages (RMS) : Vp = {self.Vp_rms_applied:.4f} V   Vs = {self.Vs_rms:.4f} V",
+            f"  Currents (RMS) : Ip = {self.Ip_rms:.4f} A   Is = {self.Is_rms:.4f} A",
+            f"  Load phase     : {math.degrees(self.load_phase_rad):.1f} deg",
+            f"  Power          : P_out={self.P_out_W:.4f} W   P_cu={self.P_cu_W:.4f} W   "
+            f"P_core={self.P_core_W:.4f} W   P_in={self.P_in_W:.4f} W",
+            f"  Efficiency     : {self.efficiency*100:.2f}%",
+        ]
+        if self.violations:
+            lines.append(f"  *** CONSTRAINT VIOLATIONS: {', '.join(self.violations)} ***")
+        else:
+            lines.append("  Constraints    : all satisfied")
+        return "\n".join(lines)
+
+
+# ---------------------------------------------------------------------------
+# Constraints dataclass
+# ---------------------------------------------------------------------------
+
+@dataclass
+class SimConstraints:
+    """
+    Operating limits.  Set any limit to float('inf') / None to disable it.
+
+    Parameters
+    ----------
+    B_max_T : float
+        Maximum peak flux density (Tesla).  Default 0.3 T.
+    Vp_max : float
+        Maximum primary RMS voltage.
+    Vs_max : float
+        Maximum secondary RMS voltage across load.
+    Ip_max : float
+        Maximum primary RMS current.
+    Is_max : float
+        Maximum secondary RMS current.
+    P_out_max_W : float
+        Maximum output power (W).
+    """
+    B_max_T: float = 0.3
+    Vp_max: float = float('inf')
+    Vs_max: float = float('inf')
+    Ip_max: float = float('inf')
+    Is_max: float = float('inf')
+    P_out_max_W: float = float('inf')
+
+
+# ---------------------------------------------------------------------------
+# Simulator
+# ---------------------------------------------------------------------------
+
+class ToroidSimulator:
+    """
+    Simulate an operating point of a toroidal transformer designed with
+    designer.py.
+
+    Parameters
+    ----------
+    core : ToroidCore
+        Core geometry and magnetic material parameters.  Magnetic fields
+        (Ae_mm2, Ve_mm3, Pv_func) are read directly from the core object;
+        if omitted they default to the geometric values.
+    primary : WindingResult
+        Output of analyse_winding / design_transformer for the primary.
+    secondary : WindingResult
+        Output of analyse_winding / design_transformer for the secondary.
+    """
+
+    def __init__(
+        self,
+        core: ToroidCore,
+        primary: WindingResult,
+        secondary: WindingResult,
+    ) -> None:
+        self.core = core
+        self.primary = primary
+        self.secondary = secondary
+
+    # ------------------------------------------------------------------
+    # Internal helpers
+    # ------------------------------------------------------------------
+
+    def _effective_turns(self, winding: WindingResult, tap: int) -> int:
+        """Cumulative turns up to and including tap (1-indexed)."""
+        n_taps = len(winding.spec.taps) - 1   # number of valid tap numbers
+        if tap < 1 or tap > n_taps:
+            raise ValueError(
+                f"Tap {tap} out of range 1..{n_taps} for winding '{winding.spec.name}'"
+            )
+        return winding.spec.cumulative_turns()[tap - 1]
+
+    def _winding_resistance(self, winding: WindingResult, tap: int) -> float:
+        """Total DC resistance for turns up to (and including) tap (ohms)."""
+        n_taps = len(winding.spec.taps) - 1
+        if tap < 1 or tap > n_taps:
+            raise ValueError(
+                f"Tap {tap} out of range 1..{n_taps} for winding '{winding.spec.name}'"
+            )
+        # Sum resistance across segments 0..tap-1
+        return sum(winding.segments[i].resistance_ohm for i in range(tap))
+
+    def _effective_Ae_mm2(self) -> float:
+        """Ae_mm2: use core field if set, else geometric cross-section."""
+        if self.core.Ae_mm2 is not None:
+            return self.core.Ae_mm2
+        return self.core.cross_section_area_mm2
+
+    def _effective_Ve_mm3(self) -> float:
+        """Ve_mm3: use core field if set, else Ae_mm2 * mean_path_length_mm."""
+        if self.core.Ve_mm3 is not None:
+            return self.core.Ve_mm3
+        return self._effective_Ae_mm2() * self.core.mean_path_length_mm
+
+    def _core_loss_W(self, B_peak_T: float, freq_hz: float) -> float:
+        """Compute core loss in watts.  Pv_func returns kW/m³; Ve in mm³."""
+        Ve_mm3 = self._effective_Ve_mm3()
+        if Ve_mm3 <= 0:
+            return 0.0
+
+        Ve_m3 = Ve_mm3 * 1e-9  # mm³ → m³
+
+        if self.core.Pv_func is not None:
+            # User-supplied function returns kW/m³
+            try:
+                Pv_kW_m3 = float(self.core.Pv_func(freq_hz, B_peak_T))
+            except Exception:
+                return 0.0
+            return Pv_kW_m3 * 1e3 * Ve_m3   # kW/m³ → W/m³, then × m³ = W
+
+        # Built-in fallback: core_loss_Pv also returns kW/m³
+        try:
+            from core import core_loss_Pv
+            Pv_kW_m3 = float(core_loss_Pv(B_peak_T, freq_hz))
+            return Pv_kW_m3 * 1e3 * Ve_m3
+        except Exception:
+            return 0.0
+
+    # ------------------------------------------------------------------
+    # Main simulate method
+    # ------------------------------------------------------------------
+
+    def simulate(
+        self,
+        Vp_rms: float,
+        freq_hz: float,
+        primary_tap: int,
+        secondary_tap: int,
+        Z_load: Tuple[float, float] = (0.0, 0.0),
+        constraints: Optional[SimConstraints] = None,
+    ) -> SimResult:
+        """
+        Compute one operating point.
+
+        Parameters
+        ----------
+        Vp_rms : float
+            Applied primary RMS voltage (V).
+        freq_hz : float
+            Excitation frequency (Hz).
+        primary_tap : int
+            Primary tap number (1-indexed).
+        secondary_tap : int
+            Secondary tap number (1-indexed).
+        Z_load : (R, X) tuple
+            Complex load impedance in ohms.  R is the resistive part,
+            X is the reactive part (positive = inductive, negative = capacitive).
+            For a purely resistive load use (R, 0.0).
+            For a complex load the simulator computes real output power as
+            P_out = |Is|² · R_load.
+        constraints : SimConstraints, optional
+            Operating limits.  If None, no constraint checking is performed.
+
+        Returns
+        -------
+        SimResult
+        """
+        if Vp_rms <= 0:
+            raise ValueError("Vp_rms must be > 0")
+        if freq_hz <= 0:
+            raise ValueError("freq_hz must be > 0")
+
+        R_load, X_load = float(Z_load[0]), float(Z_load[1])
+        Z_load_complex = complex(R_load, X_load)
+
+        if abs(Z_load_complex) == 0:
+            raise ValueError("Z_load magnitude must be > 0 (short circuit not supported)")
+
+        # --- Turns and resistance ------------------------------------------
+        Np = self._effective_turns(self.primary, primary_tap)
+        Ns = self._effective_turns(self.secondary, secondary_tap)
+        Rp = self._winding_resistance(self.primary, primary_tap)
+        Rs = self._winding_resistance(self.secondary, secondary_tap)
+
+        if Np == 0 or Ns == 0:
+            raise ValueError("Effective turns cannot be zero")
+
+        turns_ratio = Ns / Np           # a = Ns/Np
+
+        # --- Flux density ---------------------------------------------------
+        # The primary voltage applied to the winding minus the resistive drop
+        # sets the flux. We iterate once: compute ideal Ip, correct Vp seen by
+        # core, then recompute.
+        #
+        # Iteration:
+        #   Pass 0 (ideal): B from full Vp
+        #   Pass 1: B from (Vp - Ip_0 * Rp)  where Ip_0 is ideal primary current
+        #
+        Ae_m2 = self._effective_Ae_mm2() * 1e-6
+
+        def _compute_op(Vp_core: float):
+            """Compute operating point given effective core voltage."""
+            B = Vp_core / (4.44 * Np * Ae_m2 * freq_hz)
+            # Ideal open-circuit secondary voltage
+            Vs_oc = complex(Vp_core * turns_ratio, 0.0)
+            # Secondary loop: Vs_oc = Is * (Rs + Z_load)
+            Z_sec_total = complex(Rs, 0.0) + Z_load_complex
+            Is_phasor = Vs_oc / Z_sec_total
+            # Voltage across load
+            Vs_phasor = Is_phasor * Z_load_complex
+            # Reflect secondary current to primary (ideal transformer)
+            Ip_reflected = Is_phasor / turns_ratio   # phasor, Amperes
+            return B, Is_phasor, Vs_phasor, Ip_reflected
+
+        # Pass 0: ideal (ignore primary drop for now)
+        _, Is0, _, Ip0 = _compute_op(Vp_rms)
+        Ip0_mag = abs(Ip0)
+
+        # Pass 1: correct for primary winding voltage drop
+        # Primary drop is Ip * Rp (in phase with current; Rp is real)
+        # Vp_core ≈ Vp_rms - Ip0 * Rp  (phasor subtraction)
+        # For simplicity treat Ip0 as real-valued magnitude for the correction
+        # (conservative: subtracts in phase with voltage)
+        Vp_core = max(Vp_rms - Ip0_mag * Rp, 0.0)
+
+        B_peak_T, Is_phasor, Vs_phasor, Ip_phasor = _compute_op(Vp_core)
+
+        Is_rms = abs(Is_phasor)
+        Vs_rms_out = abs(Vs_phasor)
+        Ip_rms = abs(Ip_phasor)
+
+        # Load phase angle (angle of Is relative to Vs across load)
+        load_phase_rad = cmath.phase(Is_phasor) - cmath.phase(Vs_phasor)
+
+        # --- Power -----------------------------------------------------------
+        # Real power out = |Is|² · R_load
+        P_out = Is_rms ** 2 * R_load
+        P_cu_p = Ip_rms ** 2 * Rp
+        P_cu_s = Is_rms ** 2 * Rs
+        P_cu = P_cu_p + P_cu_s
+        P_core = self._core_loss_W(B_peak_T, freq_hz)
+        P_in = P_out + P_cu + P_core
+        efficiency = P_out / P_in if P_in > 0 else 1.0
+
+        # --- Constraint checking --------------------------------------------
+        violations: list[str] = []
+        if constraints is not None:
+            if B_peak_T > constraints.B_max_T:
+                violations.append(
+                    f"B_peak={B_peak_T*1e3:.1f}mT > B_max={constraints.B_max_T*1e3:.1f}mT"
+                )
+            if Vp_rms > constraints.Vp_max:
+                violations.append(
+                    f"Vp={Vp_rms:.2f}V > Vp_max={constraints.Vp_max:.2f}V"
+                )
+            if Vs_rms_out > constraints.Vs_max:
+                violations.append(
+                    f"Vs={Vs_rms_out:.2f}V > Vs_max={constraints.Vs_max:.2f}V"
+                )
+            if Ip_rms > constraints.Ip_max:
+                violations.append(
+                    f"Ip={Ip_rms:.4f}A > Ip_max={constraints.Ip_max:.4f}A"
+                )
+            if Is_rms > constraints.Is_max:
+                violations.append(
+                    f"Is={Is_rms:.4f}A > Is_max={constraints.Is_max:.4f}A"
+                )
+            if P_out > constraints.P_out_max_W:
+                violations.append(
+                    f"P_out={P_out:.2f}W > P_out_max={constraints.P_out_max_W:.2f}W"
+                )
+
+        return SimResult(
+            Vp_rms=Vp_rms,
+            freq_hz=freq_hz,
+            primary_tap=primary_tap,
+            secondary_tap=secondary_tap,
+            Np_eff=Np,
+            Ns_eff=Ns,
+            turns_ratio=turns_ratio,
+            B_peak_T=B_peak_T,
+            Vs_rms=Vs_rms_out,
+            Vp_rms_applied=Vp_rms,
+            Ip_rms=Ip_rms,
+            Is_rms=Is_rms,
+            load_phase_rad=load_phase_rad,
+            Rp=Rp,
+            Rs=Rs,
+            P_out_W=P_out,
+            P_cu_primary_W=P_cu_p,
+            P_cu_secondary_W=P_cu_s,
+            P_cu_W=P_cu,
+            P_core_W=P_core,
+            P_in_W=P_in,
+            efficiency=efficiency,
+            violations=violations,
+        )
+
+    def optimize(
+        self,
+        freq_hz: float,
+        Z_load: Tuple[float, float],
+        target_power_W: float,
+        constraints: SimConstraints,
+        Vp_min: float = 1.0,
+        Vp_max: float = 50.0,
+        Vp_steps: int = 100,
+        power_tol_pct: float = 2.0,
+    ) -> Optional["OptimizeResult"]:
+        """
+        Find the tap combination and input voltage that best delivers
+        ``target_power_W`` to the load while satisfying all constraints.
+
+        The search evaluates every (primary_tap, secondary_tap, Vp) triple.
+        A candidate is *feasible* if it passes all hard limits in
+        ``constraints`` (B_max_T, Vp_max, Vs_max, Ip_max, Is_max).
+        ``constraints.P_out_max_W`` is treated as a hard upper bound too.
+
+        Among feasible candidates:
+
+        1. If any candidate delivers power within ``power_tol_pct`` % of
+           ``target_power_W``, the one with the highest efficiency is returned.
+        2. Otherwise (fallback), the candidate with power closest to
+           ``target_power_W`` is returned; ties broken by highest efficiency.
+
+        Parameters
+        ----------
+        freq_hz : float
+            Operating frequency (Hz).
+        Z_load : (R, X) tuple
+            Complex load impedance (ohms).
+        target_power_W : float
+            Desired output power (W).
+        constraints : SimConstraints
+            Hard operating limits.  All fields are enforced.
+        Vp_min, Vp_max : float
+            Input voltage search range (V).
+        Vp_steps : int
+            Number of voltage steps (linearly spaced, inclusive of endpoints).
+        power_tol_pct : float
+            Acceptable power delivery error as % of target.  Default 2 %.
+
+        Returns
+        -------
+        OptimizeResult, or None if no feasible point exists at all.
+        """
+        import numpy as np
+
+        if target_power_W <= 0:
+            raise ValueError("target_power_W must be > 0")
+        if Vp_min >= Vp_max:
+            raise ValueError("Vp_min must be < Vp_max")
+        if Vp_steps < 2:
+            raise ValueError("Vp_steps must be >= 2")
+
+        voltages = np.linspace(Vp_min, Vp_max, Vp_steps)
+        n_ptaps = len(self.primary.spec.taps) - 1
+        n_staps = len(self.secondary.spec.taps) - 1
+
+        tol_abs = target_power_W * power_tol_pct / 100.0
+
+        # Best within tolerance
+        best_in_tol: Optional[SimResult] = None
+        best_in_tol_eff = -1.0
+
+        # Best fallback (closest power, then best efficiency)
+        best_fallback: Optional[SimResult] = None
+        best_fallback_err = float('inf')
+        best_fallback_eff = -1.0
+
+        for p_tap in range(1, n_ptaps + 1):
+            for s_tap in range(1, n_staps + 1):
+                for Vp in voltages:
+                    try:
+                        r = self.simulate(
+                            Vp_rms=float(Vp),
+                            freq_hz=freq_hz,
+                            primary_tap=p_tap,
+                            secondary_tap=s_tap,
+                            Z_load=Z_load,
+                            constraints=constraints,
+                        )
+                    except ValueError:
+                        continue
+
+                    # Skip if any hard constraint violated
+                    if not r.feasible:
+                        continue
+
+                    err = abs(r.P_out_W - target_power_W)
+
+                    if err <= tol_abs:
+                        if r.efficiency > best_in_tol_eff:
+                            best_in_tol_eff = r.efficiency
+                            best_in_tol = r
+                    else:
+                        if (err < best_fallback_err or
+                                (err == best_fallback_err and r.efficiency > best_fallback_eff)):
+                            best_fallback_err = err
+                            best_fallback_eff = r.efficiency
+                            best_fallback = r
+
+        if best_in_tol is not None:
+            return OptimizeResult(result=best_in_tol, target_power_W=target_power_W,
+                                  power_tol_pct=power_tol_pct, met_target=True)
+        if best_fallback is not None:
+            return OptimizeResult(result=best_fallback, target_power_W=target_power_W,
+                                  power_tol_pct=power_tol_pct, met_target=False)
+        return None
+
+
+# ---------------------------------------------------------------------------
+# Optimize result dataclass
+# ---------------------------------------------------------------------------
+
+@dataclass
+class OptimizeResult:
+    """
+    Result of ToroidSimulator.optimize().
+
+    Attributes
+    ----------
+    result : SimResult
+        The best operating point found.
+    target_power_W : float
+        The requested target output power.
+    power_tol_pct : float
+        The tolerance that was used (%).
+    met_target : bool
+        True if result.P_out_W is within power_tol_pct % of target_power_W.
+        False means the optimizer returned the closest feasible point it
+        could find (fallback mode).
+    """
+    result: SimResult
+    target_power_W: float
+    power_tol_pct: float
+    met_target: bool
+
+    @property
+    def power_error_pct(self) -> float:
+        return abs(self.result.P_out_W - self.target_power_W) / self.target_power_W * 100.0
+
+    def __str__(self) -> str:
+        status = "TARGET MET" if self.met_target else "FALLBACK (closest feasible)"
+        lines = [
+            f"=== Optimize result [{status}] ===",
+            f"  Target power   : {self.target_power_W:.3f} W  "
+            f"(tol={self.power_tol_pct:.1f}%,  actual error={self.power_error_pct:.2f}%)",
+            str(self.result),
+        ]
+        return "\n".join(lines)
+
+
+# ---------------------------------------------------------------------------
+# Multi-point sweep: optimise over a grid of frequencies × loads
+# ---------------------------------------------------------------------------
+
+@dataclass
+class SweepEntry:
+    """
+    One row of the operating-point sweep dataset.
+
+    Fields that come from the optimizer search inputs are prefixed with no
+    special marker.  Fields derived from the best SimResult are included
+    directly so the object is self-contained and easy to serialise.
+    """
+    # --- Search inputs -------------------------------------------------------
+    freq_hz: float
+    Z_load_R: float             # resistive part of load (ohms)
+    Z_load_X: float             # reactive part of load (ohms)
+    target_power_W: float
+
+    # --- Optimizer outcome ---------------------------------------------------
+    met_target: bool            # True = within tolerance; False = fallback
+    power_error_pct: float      # |P_out - target| / target * 100
+
+    # --- Best operating point (None fields mean no feasible solution) --------
+    primary_tap: Optional[int]
+    secondary_tap: Optional[int]
+    Vp_rms: Optional[float]
+    Np_eff: Optional[int]
+    Ns_eff: Optional[int]
+    turns_ratio: Optional[float]
+    B_peak_T: Optional[float]
+    Vs_rms: Optional[float]
+    Ip_rms: Optional[float]
+    Is_rms: Optional[float]
+    load_phase_deg: Optional[float]
+    Rp_ohm: Optional[float]
+    Rs_ohm: Optional[float]
+    P_out_W: Optional[float]
+    P_cu_W: Optional[float]
+    P_cu_primary_W: Optional[float]
+    P_cu_secondary_W: Optional[float]
+    P_core_W: Optional[float]
+    P_in_W: Optional[float]
+    efficiency: Optional[float]
+
+    # CSV column order (class-level constant)
+    _CSV_FIELDS: ClassVar[list[str]] = [
+        "freq_hz", "Z_load_R", "Z_load_X", "target_power_W",
+        "met_target", "power_error_pct",
+        "primary_tap", "secondary_tap", "Vp_rms",
+        "Np_eff", "Ns_eff", "turns_ratio",
+        "B_peak_T", "Vs_rms", "Ip_rms", "Is_rms", "load_phase_deg",
+        "Rp_ohm", "Rs_ohm",
+        "P_out_W", "P_cu_W", "P_cu_primary_W", "P_cu_secondary_W",
+        "P_core_W", "P_in_W", "efficiency",
+    ]
+
+    def as_dict(self) -> dict:
+        return {f: getattr(self, f) for f in self._CSV_FIELDS}
+
+    @staticmethod
+    def _no_solution(
+        freq_hz: float,
+        Z_load: Tuple[float, float],
+        target_power_W: float,
+    ) -> "SweepEntry":
+        """Create a row representing an infeasible / no-solution point."""
+        return SweepEntry(
+            freq_hz=freq_hz, Z_load_R=Z_load[0], Z_load_X=Z_load[1],
+            target_power_W=target_power_W,
+            met_target=False, power_error_pct=float('nan'),
+            primary_tap=None, secondary_tap=None, Vp_rms=None,
+            Np_eff=None, Ns_eff=None, turns_ratio=None,
+            B_peak_T=None, Vs_rms=None, Ip_rms=None, Is_rms=None,
+            load_phase_deg=None, Rp_ohm=None, Rs_ohm=None,
+            P_out_W=None, P_cu_W=None, P_cu_primary_W=None,
+            P_cu_secondary_W=None, P_core_W=None, P_in_W=None,
+            efficiency=None,
+        )
+
+    @staticmethod
+    def _from_opt(
+        freq_hz: float,
+        Z_load: Tuple[float, float],
+        target_power_W: float,
+        opt: "OptimizeResult",
+    ) -> "SweepEntry":
+        r = opt.result
+        return SweepEntry(
+            freq_hz=freq_hz, Z_load_R=Z_load[0], Z_load_X=Z_load[1],
+            target_power_W=target_power_W,
+            met_target=opt.met_target,
+            power_error_pct=opt.power_error_pct,
+            primary_tap=r.primary_tap,
+            secondary_tap=r.secondary_tap,
+            Vp_rms=r.Vp_rms,
+            Np_eff=r.Np_eff,
+            Ns_eff=r.Ns_eff,
+            turns_ratio=r.turns_ratio,
+            B_peak_T=r.B_peak_T,
+            Vs_rms=r.Vs_rms,
+            Ip_rms=r.Ip_rms,
+            Is_rms=r.Is_rms,
+            load_phase_deg=math.degrees(r.load_phase_rad),
+            Rp_ohm=r.Rp,
+            Rs_ohm=r.Rs,
+            P_out_W=r.P_out_W,
+            P_cu_W=r.P_cu_W,
+            P_cu_primary_W=r.P_cu_primary_W,
+            P_cu_secondary_W=r.P_cu_secondary_W,
+            P_core_W=r.P_core_W,
+            P_in_W=r.P_in_W,
+            efficiency=r.efficiency,
+        )
+
+
+def sweep_operating_points(
+    sim: ToroidSimulator,
+    frequencies: list[float],
+    loads: list[Tuple[float, float]],
+    target_power_W: float,
+    constraints: SimConstraints,
+    Vp_min: float = 1.0,
+    Vp_max: float = 50.0,
+    Vp_steps: int = 100,
+    power_tol_pct: float = 2.0,
+) -> list[SweepEntry]:
+    """
+    Optimise over a grid of frequencies × complex loads.
+
+    For each (freq, load) combination, calls ``sim.optimize()`` to find the
+    tap combination and input voltage that best delivers ``target_power_W``
+    while satisfying all constraints.
+
+    Parameters
+    ----------
+    sim : ToroidSimulator
+    frequencies : list[float]
+        Frequencies to sweep (Hz).
+    loads : list[(R, X)]
+        Complex load impedances to sweep (ohms).
+    target_power_W : float
+        Desired output power for every operating point.
+    constraints : SimConstraints
+        Hard limits applied at every point.
+    Vp_min, Vp_max : float
+        Input voltage search range (V).
+    Vp_steps : int
+        Number of voltage steps per optimize call.
+    power_tol_pct : float
+        Acceptable power delivery error (%).
+
+    Returns
+    -------
+    list[SweepEntry]
+        One entry per (freq, load) pair, in the order
+        frequencies × loads (freq varies slowest).
+    """
+    entries: list[SweepEntry] = []
+    for freq in frequencies:
+        for Z in loads:
+            opt = sim.optimize(
+                freq_hz=freq,
+                Z_load=Z,
+                target_power_W=target_power_W,
+                constraints=constraints,
+                Vp_min=Vp_min,
+                Vp_max=Vp_max,
+                Vp_steps=Vp_steps,
+                power_tol_pct=power_tol_pct,
+            )
+            if opt is None:
+                entries.append(SweepEntry._no_solution(freq, Z, target_power_W))
+            else:
+                entries.append(SweepEntry._from_opt(freq, Z, target_power_W, opt))
+    return entries
+
+
+def sweep_to_csv(entries: list[SweepEntry], path: str) -> None:
+    """
+    Write a sweep dataset to a CSV file.
+
+    Parameters
+    ----------
+    entries : list[SweepEntry]
+        Output of ``sweep_operating_points()``.
+    path : str
+        Destination file path (created or overwritten).
+    """
+    import csv
+    if not entries:
+        return
+    fields = SweepEntry._CSV_FIELDS
+    with open(path, "w", newline="") as fh:
+        writer = csv.DictWriter(fh, fieldnames=fields)
+        writer.writeheader()
+        for e in entries:
+            writer.writerow(e.as_dict())
+
+
+# ---------------------------------------------------------------------------
+# Convenience: sweep taps at fixed Vp / freq
+# ---------------------------------------------------------------------------
+
+def sweep_taps(
+    sim: ToroidSimulator,
+    Vp_rms: float,
+    freq_hz: float,
+    Z_load: Tuple[float, float],
+    constraints: Optional[SimConstraints] = None,
+) -> list[SimResult]:
+    """
+    Evaluate all primary × secondary tap combinations for a fixed operating
+    point.
+
+    Returns a list of SimResult, sorted by descending efficiency (feasible
+    results first).
+    """
+    n_ptaps = len(sim.primary.spec.taps) - 1
+    n_staps = len(sim.secondary.spec.taps) - 1
+    results = []
+    for p in range(1, n_ptaps + 1):
+        for s in range(1, n_staps + 1):
+            try:
+                r = sim.simulate(
+                    Vp_rms=Vp_rms,
+                    freq_hz=freq_hz,
+                    primary_tap=p,
+                    secondary_tap=s,
+                    Z_load=Z_load,
+                    constraints=constraints,
+                )
+                results.append(r)
+            except ValueError:
+                continue
+
+    # Sort: feasible first, then by descending efficiency
+    results.sort(key=lambda r: (not r.feasible, -r.efficiency))
+    return results
+
+
+# ---------------------------------------------------------------------------
+# Example / smoke-test
+# ---------------------------------------------------------------------------
+def pcv(f_hz, B_T):
+    pcv_ref = 25    # kW/m^3
+    a = 1.3
+    b = 2.8
+
+    ret = pcv_ref * (f_hz / 20000)**a * (B_T / 0.2)**b
+    return ret
+
+
+if __name__ == "__main__":
+
+
+    from designer import WindingSpec, design_transformer
+
+    Ae_mm2 = 142.5
+    Le_mm = 106.8
+    Ve_mm3 = Ae_mm2 * Le_mm
+
+    core = ToroidCore(ID_mm=21.5, OD_mm=46.5, height_mm=22.8, Ae_mm2=Ae_mm2, Ve_mm3=Ve_mm3, Pv_func=pcv)
+
+    primary = WindingSpec(
+        awg=[22, 22],
+        taps=[0, 25, 50],
+        name="primary",
+    )
+    secondary = WindingSpec(
+        awg=[22, 22, 22, 26],
+        taps=[0, 100, 50, 50, 50],
+        name="secondary",
+    )
+
+    pri_result, sec_result = design_transformer(core, [primary, secondary])
+
+    # Ae_mm2 / Ve_mm3 / Pv_func not set on core → simulator uses geometric defaults
+    sim = ToroidSimulator(core=core, primary=pri_result, secondary=sec_result)
+
+    constraints = SimConstraints(
+        B_max_T=1.0,
+        Vp_max=50.0,
+        Vs_max=90.0,
+        Ip_max=5.0,
+        Is_max=2.0,
+        P_out_max_W=25.0,
+    )
+
+    print("=== Single operating point (10 ohm resistive load) ===")
+    result = sim.simulate(
+        Vp_rms=12.0,
+        freq_hz=256.0,
+        primary_tap=2,
+        secondary_tap=1,
+        Z_load=(100.0, 0.0),
+        constraints=constraints,
+    )
+    print(result)
+    print()
+
+    print("=== Complex load (8 ohm + 1 mH at 50 kHz -> X ~= 314 ohm) ===")
+    X = 2 * math.pi * 50_000.0 * 1e-3
+    result2 = sim.simulate(
+        Vp_rms=24.0,
+        freq_hz=50_000.0,
+        primary_tap=2,
+        secondary_tap=1,
+        Z_load=(8.0, X),
+        constraints=constraints,
+    )
+    print(result2)
+    print()
+
+    print("=== Tap sweep ===")
+    all_results = sweep_taps(sim, 24.0, 50_000.0, (10.0, 0.0), constraints)
+    for r in all_results:
+        feasible = "OK" if r.feasible else "VIOL"
+        print(
+            f"  P{r.primary_tap}/S{r.secondary_tap}  "
+            f"Np={r.Np_eff:3d} Ns={r.Ns_eff:3d}  "
+            f"Vs={r.Vs_rms:.3f}V  Is={r.Is_rms:.4f}A  "
+            f"P_out={r.P_out_W:.3f}W  eff={r.efficiency*100:.2f}%  [{feasible}]"
+        )
+
+    print()
+    print("=== Optimize: find best taps + Vp to deliver ~10 W at 256 Hz ===")
+    opt = sim.optimize(
+        freq_hz=256.0,
+        Z_load=(10.0, 0.0),
+        target_power_W=25.0,
+        constraints=constraints,
+        Vp_min=1.0,
+        Vp_max=50.0,
+        Vp_steps=200,
+        power_tol_pct=2.0,
+    )
+    if opt is None:
+        print("  No feasible solution found.")
+    else:
+        print(opt)
+
+    print()
+    print("=== Operating-point sweep (3 freqs x 3 loads) ===")
+    freqs = [256.0, 870.0, 3140.0, 8900.0]
+    loads = [(50.0, 0.0), (100.0, 0.0), (200.0, 0.0), (600.0, 0.0)]
+    entries = sweep_operating_points(
+        sim=sim,
+        frequencies=freqs,
+        loads=loads,
+        target_power_W=25.0,
+        constraints=constraints,
+        Vp_min=1.0,
+        Vp_max=50.0,
+        Vp_steps=100,
+        power_tol_pct=2.0,
+    )
+    # Print a compact table
+    print(f"  {'freq_hz':>8}  {'R':>6}  {'X':>6}  {'met':>5}  "
+          f"{'P_out':>7}  {'err%':>5}  {'eff%':>6}  {'Vp':>6}  {'tap':>5}")
+    for e in entries:
+        met = "YES" if e.met_target else "NO"
+        if e.P_out_W is not None:
+            print(f"  {e.freq_hz:>8.1f}  {e.Z_load_R:>6.1f}  {e.Z_load_X:>6.1f}  "
+                  f"{met:>5}  {e.P_out_W:>7.3f}  {e.power_error_pct:>5.2f}  "
+                  f"{e.efficiency*100:>6.2f}  {e.Vp_rms:>6.3f}  "
+                  f"P{e.primary_tap}/S{e.secondary_tap}")
+        else:
+            print(f"  {e.freq_hz:>8.1f}  {e.Z_load_R:>6.1f}  {e.Z_load_X:>6.1f}  "
+                  f"{'NO SOL':>5}")
+
+    # Write to CSV
+    sweep_to_csv(entries, "sweep_results.csv")
+    print()
+    print("  CSV written to sweep_results.csv")

start_webapp.bat

@@ -0,0 +1,9 @@
+@echo off
+echo Starting Transformer Optimizer Web App...
+echo.
+echo Open your browser to: http://localhost:5000
+echo.
+echo Press Ctrl+C to stop the server
+echo.
+python app.py
+pause