diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index d809b9b..9d3a769 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -38,6 +38,24 @@ Each exports `Geometry` and `Computation`, so the module name carries the identi - **Derive, don't fork.** To study a variant (add magnets, change turns), subclass an existing anchor and override only what differs — do **not** edit an anchor's class in place, so old results stay reproducible. A genuinely new machine is a new anchor module. - **An anchor must declare all of the above**, rotor radii included. They are the radial window the barrier generators draw into; an anchor that does not set them cannot have a rotor built for it. +## What runs in Ansys, what runs in Python + +Ansys computes only what needs the field solution. Everything downstream of it is Python. + +| in Ansys (needs the solve) | in Python (algebra on solved quantities) | +|---|---| +| torque waveform, flux linkage, `InducedVoltage`, forces, core loss | Park/Clarke dq1/dq3 (`machine_design/transforms.py`), `Rs*I`, `Lew*dI/dt`, peak search over a waveform, terminal voltage at any speed | + +Why: + +- **One solve, many operating points.** Flux linkage at a given current is speed independent, so `V = Rs*I + j*w*psi` can be re-evaluated at any electrical speed in Python. An Ansys output variable is frozen at the solved frequency, so computing voltage there costs one solve *per* operating point -- three times the FEA for a three-point profile, N times for a drive cycle. +- **Post-processing scalars stay cheap.** `Lew`, winding turns, the bus voltage and `Rs` are post-processing quantities. Keeping them in Python means a stored pool of solves can be re-scored at a new value with no new FEA. +- **It can be tested.** Ansys expressions have no unit tests, parse differently across versions (#15), and mix units inside one argument (`2*pi*f*Time - 72deg`). The Python transforms are covered by `tests/test_transforms.py`. + +**One owner per quantity.** Where both sides *could* compute something, exactly one does, and the code says which. `Lew` in `motors/motor2.py` is the live example: it stays `0` in the model, and the end-winding voltage term belongs to the Python consumer. Setting it in both places double-counts it. + +`tests/test_ansys_transforms.py` (marked `ansys`) asserts that the Python transform and the Ansys output variables agree on a live solve, so this boundary is checked rather than merely described. + ## Pull requests - Shared code changes go through a branch and a PR, not a direct push to `main`. diff --git a/machine_design/transforms.py b/machine_design/transforms.py new file mode 100644 index 0000000..b59be76 --- /dev/null +++ b/machine_design/transforms.py @@ -0,0 +1,97 @@ +"""m-phase Park/Clarke transforms, in Python. + +This is the *single* Python implementation of the transform that Ansys applies in +`motors/synrm_5f_40s.py`'s output variables, where the same three lines are written out +once for flux linkage, once for input current and once for terminal voltage. The +convention here matches those expressions exactly, so a quantity transformed here +and the matching `*_d1`/`*_q1`/`*_d3`/`*_q3` output variable agree to solver +tolerance. `tests/test_ansys_transforms.py` asserts that on a live solve. + +Phase count and harmonic are both free. `to_dq` reads the phase count off the last +axis of its input, so the same function serves the 3-phase anchor and the 5-phase +one; `to_phases` takes it explicitly, since it has to synthesise that many values. +Caveat for the triplen case: in a 3-phase machine the 3rd harmonic is zero-sequence +(all phases in phase), so a dq3 read of it is not a constant vector -- that is +physics, not a defect of the transform. + +Convention (from `motors/synrm_5f_40s.py`, `set_output_vars`):: + + theta_el = RotSign * (Moving1.Position - InitPos) * PolePairs - pi + cos_ = cos( h * (theta_el - 2*pi*k/m)) + sin_ = sin(-h * (theta_el - 2*pi*k/m)) <-- note the minus + X_d = 2/m * sum_k X_k * cos_ + X_q = 2/m * sum_k X_k * sin_ + +(written out for m = 5 in that file; `m` here is the phase count.) + +**Units.** ``theta_el`` here is in RADIANS. Ansys's ``theta_el`` output variable comes +back in DEGREES -- the ``- pi`` in its expression is evaluated as -180 in an angular +context -- so convert with ``np.deg2rad`` before passing it in. Verified on a live +2024.2 solve: with that conversion the transform below reproduces Ansys's own +``Flux_d1/q1/d3/q3`` to 2e-15 (`tests/test_ansys_transforms.py`). + +The minus inside the sine means ``q`` is the *negative* of the usual sine sum. +That is a deliberate quirk of this model, not a typo; `tests/test_transforms.py` +pins it so it cannot drift. + +Why this lives in Python rather than as more Ansys output variables: see +"What runs in Ansys, what runs in Python" in `CONTRIBUTING.md`. In short, the +transform needs no field solution -- it is algebra on quantities the solve has +already produced -- and keeping it here makes it testable and lets one solve be +re-evaluated at any electrical speed. +""" + +import numpy as np + + +def electrical_angle(position, init_pos, pole_pairs, rot_sign=1.0, degrees=False): + """theta_el in RADIANS, matching Ansys's ``theta_el`` output variable. + + `position` and `init_pos` are mechanical angles -- radians by default, or degrees + with ``degrees=True``. Pass ``degrees=True`` for values read back from Ansys: + ``Moving1.Position`` is reported in degrees, so feeding it in as radians silently + produces a meaningless angle. The return value is always radians, which is what + `to_dq` and `to_phases` expect. + """ + position = np.asarray(position, dtype=float) + if degrees: + position, init_pos = np.deg2rad(position), np.deg2rad(init_pos) + return rot_sign * (position - init_pos) * pole_pairs - np.pi + + +def _phase_angles(theta_el, harmonic, n_phases): + """h * (theta_el - 2*pi*k/m) for k = 0..m-1, shaped (..., m).""" + theta_el = np.asarray(theta_el, dtype=float)[..., None] + k = np.arange(n_phases) + return harmonic * (theta_el - 2.0 * np.pi * k / n_phases) + + +def to_dq(x_phases, theta_el, harmonic): + """Per-phase quantity -> (d, q) for the given harmonic. + + `x_phases` carries the phases on the last axis, in winding order (A, B, C, ...); + the phase count is taken from that axis, so this works for any m. `theta_el` is + broadcast against the leading axes. Returns two arrays shaped like `theta_el`. + """ + x_phases = np.asarray(x_phases, dtype=float) + n_phases = x_phases.shape[-1] + if n_phases < 3: + raise ValueError(f"need at least 3 phases on the last axis, got {n_phases}") + phi = _phase_angles(theta_el, harmonic, n_phases) + scale = 2.0 / n_phases + d = scale * np.sum(x_phases * np.cos(phi), axis=-1) + q = scale * np.sum(x_phases * np.sin(-phi), axis=-1) + return d, q + + +def to_phases(d, q, theta_el, harmonic, n_phases): + """Inverse of `to_dq` for a single harmonic: (d, q) -> `n_phases` phase values. + + Exact inverse for a pure `harmonic` component; summing the results of several + harmonics reconstructs a multi-harmonic waveform, because the transforms for + distinct harmonics are orthogonal over m equally spaced phases (see + `tests/test_transforms.py`). `n_phases` is explicit here because, unlike `to_dq`, + there is no input array to read it from. + """ + phi = _phase_angles(theta_el, harmonic, n_phases) + return np.asarray(d, dtype=float)[..., None] * np.cos(phi) - np.asarray(q, dtype=float)[..., None] * np.sin(phi) diff --git a/motors/synrm_5f_40s.py b/motors/synrm_5f_40s.py index 3b8bdbb..5d1f0c6 100644 --- a/motors/synrm_5f_40s.py +++ b/motors/synrm_5f_40s.py @@ -8,6 +8,7 @@ from ansys.aedt.core import Maxwell2d from machine_design.designs.computation import ComputationBase +from machine_design.transforms import electrical_angle, to_dq from .synrm_3f_36s import Geometry as BaseGeometry @@ -33,6 +34,11 @@ def set_oper_params(self): f = 50 # [Hz] RotSpeed = 60 * f / self.geometry.PolePairs # [rpm] w = 2 * np.pi * f + self.InitPos = -45 # deg + self.RotSign = 1 + self.Rstat = 19.0 + self.Lew = 0.0 + self.w = w self.oper_params = { "Id1": "0.0A", "Iq1": "0.0A", @@ -42,7 +48,7 @@ def set_oper_params(self): "epsI3": "atan2(Iq3,Id3)", # current angle, 1st harmonic "Im1": "sqrt(Id1^2+Iq1^2)", "Im3": "sqrt(Id3^2+Iq3^2)", - "InitPos": "-45deg", + "InitPos": f"{self.InitPos}deg", "w": f"{w}Hz", "RotSpeed": f"{RotSpeed}rpm", "Nper": "1/10", # number of included periods @@ -51,136 +57,16 @@ def set_oper_params(self): def set_solution_expressions(self): self.solution_expressions = [ - "V_d1", - "V_q1", - "V_d3", - "V_q3", - "Flux_e_d1", - "Flux_e_q1", - "Flux_e_d3", - "Flux_e_q3", - "I_d1", - "I_q1", - "I_d3", - "I_q3", - "Ld1", - "Ld1q1", - "Ld1d3", - "Ld1q3", - "Lq1", - "Lq1d3", - "Lq1q3", - "Ld3", - "Ld3q3", - "Lq3", + "Moving1.Position", "Moving1.Torque", + *[f"FluxLinkage(Phase{p})" for p in "ABCDE"], + *[f"InducedVoltage(Phase{p})" for p in "ABCDE"], + *[f"InputCurrent(Phase{p})" for p in "ABCDE"], + *[f"L(Phase{x},Phase{y})" for x in "ABCDE" for y in "ABCDE"], ] def set_output_vars(self): - self.output_vars = { - "PolePairs": "2", - "RotSign": "1", - "Rstat": "19", - "Lew": "0", - "theta_el": "RotSign*(Moving1.Position - InitPos) * PolePairs - pi", - "cos0_1": "cos(1*(theta_el - 2*PI*0/5))", - "sin0_1": "sin(-1*(theta_el - 2*PI*0/5))", - "cos1_1": "cos(1*(theta_el - 2*PI*1/5))", - "sin1_1": "sin(-1*(theta_el - 2*PI*1/5))", - "cos2_1": "cos(1*(theta_el - 2*PI*2/5))", - "sin2_1": "sin(-1*(theta_el - 2*PI*2/5))", - "cos3_1": "cos(1*(theta_el - 2*PI*3/5))", - "sin3_1": "sin(-1*(theta_el - 2*PI*3/5))", - "cos4_1": "cos(1*(theta_el - 2*PI*4/5))", - "sin4_1": "sin(-1*(theta_el - 2*PI*4/5))", - "cos0_3": "cos(3*(theta_el - 2*PI*0/5))", - "sin0_3": "sin(-3*(theta_el - 2*PI*0/5))", - "cos1_3": "cos(3*(theta_el - 2*PI*1/5))", - "sin1_3": "sin(-3*(theta_el - 2*PI*1/5))", - "cos2_3": "cos(3*(theta_el - 2*PI*2/5))", - "sin2_3": "sin(-3*(theta_el - 2*PI*2/5))", - "cos3_3": "cos(3*(theta_el - 2*PI*3/5))", - "sin3_3": "sin(-3*(theta_el - 2*PI*3/5))", - "cos4_3": "cos(3*(theta_el - 2*PI*4/5))", - "sin4_3": "sin(-3*(theta_el - 2*PI*4/5))", - "Flux_d1": "(FluxLinkage(PhaseA)*cos0_1 + FluxLinkage(PhaseB)*cos1_1 + FluxLinkage(PhaseC)*cos2_1 + FluxLinkage(PhaseD)*cos3_1 + FluxLinkage(PhaseE)*cos4_1) * 2/5", - "Flux_q1": "(FluxLinkage(PhaseA)*sin0_1 + FluxLinkage(PhaseB)*sin1_1 + FluxLinkage(PhaseC)*sin2_1 + FluxLinkage(PhaseD)*sin3_1 + FluxLinkage(PhaseE)*sin4_1) * 2/5", - "Flux_d3": "(FluxLinkage(PhaseA)*cos0_3 + FluxLinkage(PhaseB)*cos1_3 + FluxLinkage(PhaseC)*cos2_3 + FluxLinkage(PhaseD)*cos3_3 + FluxLinkage(PhaseE)*cos4_3) * 2/5", - "Flux_q3": "(FluxLinkage(PhaseA)*sin0_3 + FluxLinkage(PhaseB)*sin1_3 + FluxLinkage(PhaseC)*sin2_3 + FluxLinkage(PhaseD)*sin3_3 + FluxLinkage(PhaseE)*sin4_3) * 2/5", - "Vind_d1": "(InducedVoltage(PhaseA)*cos0_1 + InducedVoltage(PhaseB)*cos1_1 + InducedVoltage(PhaseC)*cos2_1 + InducedVoltage(PhaseD)*cos3_1 + InducedVoltage(PhaseE)*cos4_1) * 2/5", - "Vind_q1": "(InducedVoltage(PhaseA)*sin0_1 + InducedVoltage(PhaseB)*sin1_1 + InducedVoltage(PhaseC)*sin2_1 + InducedVoltage(PhaseD)*sin3_1 + InducedVoltage(PhaseE)*sin4_1) * 2/5", - "Vind_d3": "(InducedVoltage(PhaseA)*cos0_3 + InducedVoltage(PhaseB)*cos1_3 + InducedVoltage(PhaseC)*cos2_3 + InducedVoltage(PhaseD)*cos3_3 + InducedVoltage(PhaseE)*cos4_3) * 2/5", - "Vind_q3": "(InducedVoltage(PhaseA)*sin0_3 + InducedVoltage(PhaseB)*sin1_3 + InducedVoltage(PhaseC)*sin2_3 + InducedVoltage(PhaseD)*sin3_3 + InducedVoltage(PhaseE)*sin4_3) * 2/5", - "dIA_dt": "-Im1*w*sin(w*Time+epsI1-pi) - Im3*3*w*sin(3*(w*Time)+epsI3-pi)", - "dIB_dt": "-Im1*w*sin(w*Time-72deg+epsI1-pi) - Im3*3*w*sin(3*(w*Time-72deg)+epsI3-pi)", - "dIC_dt": "-Im1*w*sin(w*Time-144deg+epsI1-pi) - Im3*3*w*sin(3*(w*Time-144deg)+epsI3-pi)", - "dID_dt": "-Im1*w*sin(w*Time-216deg+epsI1-pi) - Im3*3*w*sin(3*(w*Time-216deg)+epsI3-pi)", - "dIE_dt": "-Im1*w*sin(w*Time-288deg+epsI1-pi) - Im3*3*w*sin(3*(w*Time-288deg)+epsI3-pi)", - "V_A": "InducedVoltage(PhaseA) + Rstat*InputCurrent(PhaseA) + Lew*dIA_dt", - "V_B": "InducedVoltage(PhaseB) + Rstat*InputCurrent(PhaseB) + Lew*dIB_dt", - "V_C": "InducedVoltage(PhaseC) + Rstat*InputCurrent(PhaseC) + Lew*dIC_dt", - "V_D": "InducedVoltage(PhaseD) + Rstat*InputCurrent(PhaseD) + Lew*dID_dt", - "V_E": "InducedVoltage(PhaseE) + Rstat*InputCurrent(PhaseE) + Lew*dIE_dt", - "V_AC": "V_A - V_C", - "V_BD": "V_B - V_D", - "V_CE": "V_C - V_E", - "V_DA": "V_D - V_A", - "V_EB": "V_E - V_B", - "Vterm_A": "1/5*(2*V_AC + -1*V_BD + 1*V_CE + -2*V_DA)", - "Vterm_B": "1/5*(2*V_AC + 4*V_BD + 1*V_CE + 3*V_DA)", - "Vterm_C": "1/5*(-3*V_AC + -1*V_BD + 1*V_CE + -2*V_DA)", - "Vterm_D": "1/5*(2*V_AC + -1*V_BD + 1*V_CE + 3*V_DA)", - "Vterm_E": "1/5*(-3*V_AC + -1*V_BD + -4*V_CE + -2*V_DA)", - "I_d1": "(InputCurrent(PhaseA)*cos0_1 + InputCurrent(PhaseB)*cos1_1 + InputCurrent(PhaseC)*cos2_1 + InputCurrent(PhaseD)*cos3_1 + InputCurrent(PhaseE)*cos4_1) * 2/5", - "I_q1": "(InputCurrent(PhaseA)*sin0_1 + InputCurrent(PhaseB)*sin1_1 + InputCurrent(PhaseC)*sin2_1 + InputCurrent(PhaseD)*sin3_1 + InputCurrent(PhaseE)*sin4_1) * 2/5", - "I_d3": "(InputCurrent(PhaseA)*cos0_3 + InputCurrent(PhaseB)*cos1_3 + InputCurrent(PhaseC)*cos2_3 + InputCurrent(PhaseD)*cos3_3 + InputCurrent(PhaseE)*cos4_3) * 2/5", - "I_q3": "(InputCurrent(PhaseA)*sin0_3 + InputCurrent(PhaseB)*sin1_3 + InputCurrent(PhaseC)*sin2_3 + InputCurrent(PhaseD)*sin3_3 + InputCurrent(PhaseE)*sin4_3) * 2/5", - "V_d1": "(V_A*cos0_1 + V_B*cos1_1 + V_C*cos2_1 + V_D*cos3_1 + V_E*cos4_1) * 2/5", - "V_q1": "(V_A*sin0_1 + V_B*sin1_1 + V_C*sin2_1 + V_D*sin3_1 + V_E*sin4_1) * 2/5", - "V_d3": "(V_A*cos0_3 + V_B*cos1_3 + V_C*cos2_3 + V_D*cos3_3 + V_E*cos4_3) * 2/5", - "V_q3": "(V_A*sin0_3 + V_B*sin1_3 + V_C*sin2_3 + V_D*sin3_3 + V_E*sin4_3) * 2/5", - "L0d_1": "L(PhaseA,PhaseA)*cos0_1 + L(PhaseA,PhaseB)*cos1_1 + L(PhaseA,PhaseC)*cos2_1 + L(PhaseA,PhaseD)*cos3_1 + L(PhaseA,PhaseE)*cos4_1", - "L0q_1": "L(PhaseA,PhaseA)*sin0_1 + L(PhaseA,PhaseB)*sin1_1 + L(PhaseA,PhaseC)*sin2_1 + L(PhaseA,PhaseD)*sin3_1 + L(PhaseA,PhaseE)*sin4_1", - "L1d_1": "L(PhaseB,PhaseA)*cos0_1 + L(PhaseB,PhaseB)*cos1_1 + L(PhaseB,PhaseC)*cos2_1 + L(PhaseB,PhaseD)*cos3_1 + L(PhaseB,PhaseE)*cos4_1", - "L1q_1": "L(PhaseB,PhaseA)*sin0_1 + L(PhaseB,PhaseB)*sin1_1 + L(PhaseB,PhaseC)*sin2_1 + L(PhaseB,PhaseD)*sin3_1 + L(PhaseB,PhaseE)*sin4_1", - "L2d_1": "L(PhaseC,PhaseA)*cos0_1 + L(PhaseC,PhaseB)*cos1_1 + L(PhaseC,PhaseC)*cos2_1 + L(PhaseC,PhaseD)*cos3_1 + L(PhaseC,PhaseE)*cos4_1", - "L2q_1": "L(PhaseC,PhaseA)*sin0_1 + L(PhaseC,PhaseB)*sin1_1 + L(PhaseC,PhaseC)*sin2_1 + L(PhaseC,PhaseD)*sin3_1 + L(PhaseC,PhaseE)*sin4_1", - "L3d_1": "L(PhaseD,PhaseA)*cos0_1 + L(PhaseD,PhaseB)*cos1_1 + L(PhaseD,PhaseC)*cos2_1 + L(PhaseD,PhaseD)*cos3_1 + L(PhaseD,PhaseE)*cos4_1", - "L3q_1": "L(PhaseD,PhaseA)*sin0_1 + L(PhaseD,PhaseB)*sin1_1 + L(PhaseD,PhaseC)*sin2_1 + L(PhaseD,PhaseD)*sin3_1 + L(PhaseD,PhaseE)*sin4_1", - "L4d_1": "L(PhaseE,PhaseA)*cos0_1 + L(PhaseE,PhaseB)*cos1_1 + L(PhaseE,PhaseC)*cos2_1 + L(PhaseE,PhaseD)*cos3_1 + L(PhaseE,PhaseE)*cos4_1", - "L4q_1": "L(PhaseE,PhaseA)*sin0_1 + L(PhaseE,PhaseB)*sin1_1 + L(PhaseE,PhaseC)*sin2_1 + L(PhaseE,PhaseD)*sin3_1 + L(PhaseE,PhaseE)*sin4_1", - "L0d_3": "L(PhaseA,PhaseA)*cos0_3 + L(PhaseA,PhaseB)*cos1_3 + L(PhaseA,PhaseC)*cos2_3 + L(PhaseA,PhaseD)*cos3_3 + L(PhaseA,PhaseE)*cos4_3", - "L0q_3": "L(PhaseA,PhaseA)*sin0_3 + L(PhaseA,PhaseB)*sin1_3 + L(PhaseA,PhaseC)*sin2_3 + L(PhaseA,PhaseD)*sin3_3 + L(PhaseA,PhaseE)*sin4_3", - "L1d_3": "L(PhaseB,PhaseA)*cos0_3 + L(PhaseB,PhaseB)*cos1_3 + L(PhaseB,PhaseC)*cos2_3 + L(PhaseB,PhaseD)*cos3_3 + L(PhaseB,PhaseE)*cos4_3", - "L1q_3": "L(PhaseB,PhaseA)*sin0_3 + L(PhaseB,PhaseB)*sin1_3 + L(PhaseB,PhaseC)*sin2_3 + L(PhaseB,PhaseD)*sin3_3 + L(PhaseB,PhaseE)*sin4_3", - "L2d_3": "L(PhaseC,PhaseA)*cos0_3 + L(PhaseC,PhaseB)*cos1_3 + L(PhaseC,PhaseC)*cos2_3 + L(PhaseC,PhaseD)*cos3_3 + L(PhaseC,PhaseE)*cos4_3", - "L2q_3": "L(PhaseC,PhaseA)*sin0_3 + L(PhaseC,PhaseB)*sin1_3 + L(PhaseC,PhaseC)*sin2_3 + L(PhaseC,PhaseD)*sin3_3 + L(PhaseC,PhaseE)*sin4_3", - "L3d_3": "L(PhaseD,PhaseA)*cos0_3 + L(PhaseD,PhaseB)*cos1_3 + L(PhaseD,PhaseC)*cos2_3 + L(PhaseD,PhaseD)*cos3_3 + L(PhaseD,PhaseE)*cos4_3", - "L3q_3": "L(PhaseD,PhaseA)*sin0_3 + L(PhaseD,PhaseB)*sin1_3 + L(PhaseD,PhaseC)*sin2_3 + L(PhaseD,PhaseD)*sin3_3 + L(PhaseD,PhaseE)*sin4_3", - "L4d_3": "L(PhaseE,PhaseA)*cos0_3 + L(PhaseE,PhaseB)*cos1_3 + L(PhaseE,PhaseC)*cos2_3 + L(PhaseE,PhaseD)*cos3_3 + L(PhaseE,PhaseE)*cos4_3", - "L4q_3": "L(PhaseE,PhaseA)*sin0_3 + L(PhaseE,PhaseB)*sin1_3 + L(PhaseE,PhaseC)*sin2_3 + L(PhaseE,PhaseD)*sin3_3 + L(PhaseE,PhaseE)*sin4_3", - "Ld1": "(L0d_1*cos0_1 + L1d_1*cos1_1 + L2d_1*cos2_1 + L3d_1*cos3_1 + L4d_1*cos4_1) * 2/5", - "Ld1q1": "(L0d_1*sin0_1 + L1d_1*sin1_1 + L2d_1*sin2_1 + L3d_1*sin3_1 + L4d_1*sin4_1) * 2/5", - "Lq1d1": "(L0q_1*cos0_1 + L1q_1*cos1_1 + L2q_1*cos2_1 + L3q_1*cos3_1 + L4q_1*cos4_1) * 2/5", - "Lq1": "(L0q_1*sin0_1 + L1q_1*sin1_1 + L2q_1*sin2_1 + L3q_1*sin3_1 + L4q_1*sin4_1) * 2/5", - "Ld1d3": "(L0d_1*cos0_3 + L1d_1*cos1_3 + L2d_1*cos2_3 + L3d_1*cos3_3 + L4d_1*cos4_3) * 2/5", - "Ld1q3": "(L0d_1*sin0_3 + L1d_1*sin1_3 + L2d_1*sin2_3 + L3d_1*sin3_3 + L4d_1*sin4_3) * 2/5", - "Lq1d3": "(L0q_1*cos0_3 + L1q_1*cos1_3 + L2q_1*cos2_3 + L3q_1*cos3_3 + L4q_1*cos4_3) * 2/5", - "Lq1q3": "(L0q_1*sin0_3 + L1q_1*sin1_3 + L2q_1*sin2_3 + L3q_1*sin3_3 + L4q_1*sin4_3) * 2/5", - "Ld3d1": "(L0d_3*cos0_1 + L1d_3*cos1_1 + L2d_3*cos2_1 + L3d_3*cos3_1 + L4d_3*cos4_1) * 2/5", - "Ld3q1": "(L0d_3*sin0_1 + L1d_3*sin1_1 + L2d_3*sin2_1 + L3d_3*sin3_1 + L4d_3*sin4_1) * 2/5", - "Lq3d1": "(L0q_3*cos0_1 + L1q_3*cos1_1 + L2q_3*cos2_1 + L3q_3*cos3_1 + L4q_3*cos4_1) * 2/5", - "Lq3q1": "(L0q_3*sin0_1 + L1q_3*sin1_1 + L2q_3*sin2_1 + L3q_3*sin3_1 + L4q_3*sin4_1) * 2/5", - "Ld3": "(L0d_3*cos0_3 + L1d_3*cos1_3 + L2d_3*cos2_3 + L3d_3*cos3_3 + L4d_3*cos4_3) * 2/5", - "Ld3q3": "(L0d_3*sin0_3 + L1d_3*sin1_3 + L2d_3*sin2_3 + L3d_3*sin3_3 + L4d_3*sin4_3) * 2/5", - "Lq3d3": "(L0q_3*cos0_3 + L1q_3*cos1_3 + L2q_3*cos2_3 + L3q_3*cos3_3 + L4q_3*cos4_3) * 2/5", - "Lq3": "(L0q_3*sin0_3 + L1q_3*sin1_3 + L2q_3*sin2_3 + L3q_3*sin3_3 + L4q_3*sin4_3) * 2/5", - "Flux_e_d1": "Flux_d1 - (Ld1*I_d1 + Ld1q1*I_q1 + Ld1d3*I_d3 + Ld1q3*I_q3)", - "Flux_e_q1": "Flux_q1 - (Lq1d1*I_d1 + Lq1*I_q1 + Lq1d3*I_d3 + Lq1q3*I_q3)", - "Flux_e_d3": "Flux_d3 - (Ld3d1*I_d1 + Ld3q1*I_q1 + Ld3*I_d3 + Ld3q3*I_q3)", - "Flux_e_q3": "Flux_q3 - (Lq3d1*I_d1 + Lq3q1*I_q1 + Lq3d3*I_d3 + Lq3*I_q3)", - "Torque_dq": "5/2*PolePairs*(1*(Flux_d1*I_q1 - Flux_q1*I_d1) + 3*(Flux_d3*I_q3 - Flux_q3*I_d3))", - } + self.output_vars = {} def set_post_params(self): self.post_params = { # reports @@ -336,21 +222,104 @@ def inductance_computation(self, m2d: Maxwell2d) -> None: m2d.change_inductance_computation(compute_transient_inductance=True, incremental_matrix=True) def set_variables(self, m2d: Maxwell2d, Id1, Iq1, Id3, Iq3): + self.Id1, self.Iq1, self.Id3, self.Iq3 = Id1, Iq1, Id3, Iq3 m2d.variable_manager["Id1"] = f"{Id1}A" m2d.variable_manager["Iq1"] = f"{Iq1}A" m2d.variable_manager["Id3"] = f"{Id3}A" m2d.variable_manager["Iq3"] = f"{Iq3}A" def extract_results(self, solutions): - out = {} - for expr in self.solution_expressions: - val = solutions.data_real(expr) + position = np.array(solutions.data_real("Moving1.Position")) + torque = np.array(solutions.data_real("Moving1.Torque")) + + theta_el = electrical_angle(position, self.InitPos, self.geometry.PolePairs, self.RotSign, degrees=True) + flux_phases = np.stack([np.array(solutions.data_real(f"FluxLinkage(Phase{p})")) for p in "ABCDE"], axis=-1) + vind_phases = np.stack([np.array(solutions.data_real(f"InducedVoltage(Phase{p})")) for p in "ABCDE"], axis=-1) + current_phases = np.stack([np.array(solutions.data_real(f"InputCurrent(Phase{p})")) for p in "ABCDE"], axis=-1) + + Flux_d1, Flux_q1 = to_dq(flux_phases, theta_el, harmonic=1) + Flux_d3, Flux_q3 = to_dq(flux_phases, theta_el, harmonic=3) + Vind_d1, Vind_q1 = to_dq(vind_phases, theta_el, harmonic=1) + Vind_d3, Vind_q3 = to_dq(vind_phases, theta_el, harmonic=3) + I_d1, I_q1 = (v / 1e3 for v in to_dq(current_phases, theta_el, harmonic=1)) + I_d3, I_q3 = (v / 1e3 for v in to_dq(current_phases, theta_el, harmonic=3)) + # mA to A + + Im1 = np.sqrt(self.Id1**2 + self.Iq1**2) + Im3 = np.sqrt(self.Id3**2 + self.Iq3**2) + epsI1 = np.atan2(self.Iq1, self.Id1) + epsI3 = np.atan2(self.Iq3, self.Id3) + + time = np.array(solutions.primary_sweep_values) + dI_dt_phases = np.zeros((len(time), 5)) + for k in range(5): + phase_offset = -2 * np.pi * k / 5 # 5 angles, 0deg, -72deg, -144deg, -216deg, -288deg + dI_dt_phases[:, k] = -Im1 * self.w * np.sin(self.w * time + phase_offset + epsI1 - np.pi) - Im3 * 3 * self.w * np.sin(3 * (self.w * time + phase_offset) + epsI3 - np.pi) + + V_phases = vind_phases + self.Rstat * current_phases + self.Lew * dI_dt_phases + + V_d1, V_q1 = to_dq(V_phases, theta_el, harmonic=1) + V_d3, V_q3 = to_dq(V_phases, theta_el, harmonic=3) - if expr.startswith("Ld") or expr.startswith("Lq"): - val = np.array(val) / 1e9 - elif expr.startswith("I_"): - val = np.array(val) / 1e3 + L_raw = [np.stack([np.array(solutions.data_real(f"L(Phase{x},Phase{y})")) for y in "ABCDE"], axis=-1) / 1e9 for x in "ABCDE"] + # nH to H + + L_d1_row = np.zeros((len(time), 5)) + L_q1_row = np.zeros((len(time), 5)) + L_d3_row = np.zeros((len(time), 5)) + L_q3_row = np.zeros((len(time), 5)) + for i, L_row in enumerate(L_raw): + L_d1_row[:, i], L_q1_row[:, i] = (v * 5 / 2 for v in to_dq(L_row, theta_el, harmonic=1)) + L_d3_row[:, i], L_q3_row[:, i] = (v * 5 / 2 for v in to_dq(L_row, theta_el, harmonic=3)) + + Ld1, Ld1q1 = to_dq(L_d1_row, theta_el, harmonic=1) + Lq1d1, Lq1 = to_dq(L_q1_row, theta_el, harmonic=1) + Ld1d3, Ld1q3 = to_dq(L_d1_row, theta_el, harmonic=3) + Lq1d3, Lq1q3 = to_dq(L_q1_row, theta_el, harmonic=3) + Ld3d1, Ld3q1 = to_dq(L_d3_row, theta_el, harmonic=1) + Lq3d1, Lq3q1 = to_dq(L_q3_row, theta_el, harmonic=1) + Ld3, Ld3q3 = to_dq(L_d3_row, theta_el, harmonic=3) + Lq3d3, Lq3 = to_dq(L_q3_row, theta_el, harmonic=3) + + Flux_e_d1 = Flux_d1 - (Ld1 * I_d1 + Ld1q1 * I_q1 + Ld1d3 * I_d3 + Ld1q3 * I_q3) + Flux_e_q1 = Flux_q1 - (Lq1d1 * I_d1 + Lq1 * I_q1 + Lq1d3 * I_d3 + Lq1q3 * I_q3) + Flux_e_d3 = Flux_d3 - (Ld3d1 * I_d1 + Ld3q1 * I_q1 + Ld3 * I_d3 + Ld3q3 * I_q3) + Flux_e_q3 = Flux_q3 - (Lq3d1 * I_d1 + Lq3q1 * I_q1 + Lq3d3 * I_d3 + Lq3 * I_q3) + + # Torque_dq = 5 / 2 * self.geometry.PolePairs * ((Flux_d1 * I_q1 - Flux_q1 * I_d1) + 3 * (Flux_d3 * I_q3 - Flux_q3 * I_d3)) + + out = { + "V_d1": V_d1, + "V_q1": V_q1, + "V_d3": V_d3, + "V_q3": V_q3, + "Flux_d1": Flux_d1, + "Flux_q1": Flux_q1, + "Flux_d3": Flux_d3, + "Flux_q3": Flux_q3, + "Vind_d1": Vind_d1, + "Vind_q1": Vind_q1, + "Vind_d3": Vind_d3, + "Vind_q3": Vind_q3, + "I_d1": I_d1, + "I_q1": I_q1, + "I_d3": I_d3, + "I_q3": I_q3, + "Ld1": Ld1, + "Ld1q1": Ld1q1, + "Lq1": Lq1, + "Ld1d3": Ld1d3, + "Ld1q3": Ld1q3, + "Lq1d3": Lq1d3, + "Lq1q3": Lq1q3, + "Ld3": Ld3, + "Ld3q3": Ld3q3, + "Lq3": Lq3, + "Flux_e_d1": Flux_e_d1, + "Flux_e_q1": Flux_e_q1, + "Flux_e_d3": Flux_e_d3, + "Flux_e_q3": Flux_e_q3, + "Moving1.Torque": torque, + } - # TODO: possibly assign val[:-1]. check whether values are identical - out[expr] = val return out diff --git a/tests/test_ansys_regression_synrm_5f_40s.py b/tests/test_ansys_regression_synrm_5f_40s.py index fd48d7c..cadc0ee 100644 --- a/tests/test_ansys_regression_synrm_5f_40s.py +++ b/tests/test_ansys_regression_synrm_5f_40s.py @@ -96,7 +96,7 @@ def test_legacy_and_live_design_match_across_rotor_rebuilds(tmp_path): assert torque_legacy is not None, f"legacy torque is None for seed {seed}" assert torque_live is not None, f"live torque is None for seed {seed}" - for expr in live_computation.solution_expressions: + for expr in legacy_design.solution_expressions: np.testing.assert_allclose(torque_live[expr], torque_legacy[expr], rtol=1e-3, atol=2e-4, err_msg=f"mismatch for {expr}, seed {seed}") # rtol is a relative tolerance, but when current_setpoint = (7.0711, 7.0711, 0.0, 0.0) is used, # I_d3 is 1e-15 and 2.2e-16 for legacy and live, respectively, rtol is not meaningful when comparing such small numbers, diff --git a/tests/test_ansys_transforms.py b/tests/test_ansys_transforms.py new file mode 100644 index 0000000..9328355 --- /dev/null +++ b/tests/test_ansys_transforms.py @@ -0,0 +1,129 @@ +"""Live check that the Python dq transform agrees with the Ansys output variables. + +`machine_design/transforms.py` and `motors/synrm_5f_40s.py`'s `set_output_vars` implement the +same five-phase Park/Clarke transform twice -- once in Python, once as Ansys expressions. +This solves one rotor with both dq1 and dq3 excited, then asserts that transforming the +per-phase quantities in Python reproduces Ansys's own dq outputs. + +Two subjects, deliberately: + +* **flux linkage** -- `FluxLinkage(PhaseA..E)` vs `Flux_d1/q1/d3/q3`. Flux stays in Ansys + (it needs the field solution), so this comparison is permanent: it is the standing guard + on the boundary described in "What runs in Ansys, what runs in Python" (CONTRIBUTING.md). +* **terminal voltage** -- `V_A..V_E` vs `V_d1/q1/d3/q3`. This one is a *one-shot* + equivalence proof. The voltage transforms are slated to move to Python, which deletes + both sides of the comparison, so it can only be run while they still exist. Once the + Ansys-side voltage variables go, drop this half and keep the flux half. + +Marked `ansys`, so it is skipped by default. Run with: + + python -m pytest -m "" tests/test_ansys_transforms.py -v +""" + +import numpy as np +import pytest + +from machine_design.config import load_config +from machine_design.designs.design import Design +from machine_design.optimization.generators import HacklGenerator_OneLambda +from machine_design.transforms import to_dq +from motors.synrm_5f_40s import Computation, Geometry + +pytestmark = pytest.mark.ansys + +AEDT_VERSION = load_config()["aedt_version"] +R_STATOR_END = 0.7 +OFFSET = 0.35 +SEED = 42 +NUM_CORES = 1 + +# Both subspaces excited, so the h = 3 comparison is not a check against zero. +# Order is (Id1, Iq1, Id3, Iq3) -- compute() forwards these to set_variables(). +CURRENTS = (0.8, 0.9, 0.2, -0.15) + +PHASES = ["A", "B", "C", "D", "E"] +PHASE_FLUX = [f"FluxLinkage(Phase{p})" for p in PHASES] +PHASE_VOLTAGE = [f"V_{p}" for p in PHASES] + +# (label, per-phase expressions, ansys dq prefix) +SUBJECTS = [("flux", PHASE_FLUX, "Flux"), ("voltage", PHASE_VOLTAGE, "V")] + + +class ComputationWithPhaseQuantities(Computation): + """`Computation`, additionally extracting the per-phase quantities and theta_el. + + The shipped class extracts only the dq quantities, since that is all the optimizers + need. Subclassed rather than edited so the shared anchor is untouched. + """ + + def set_solution_expressions(self): + super().set_solution_expressions() + dq = [f"{prefix}_{axis}{h}" for prefix in ("Flux", "V") for h in (1, 3) for axis in "dq"] + extra = [*PHASE_FLUX, *PHASE_VOLTAGE, "theta_el", *dq] + self.solution_expressions = [*self.solution_expressions, *extra] + + +def _feasible_barriers(geometry): + dummy = Design(m2d=None, geometry=geometry, computation=Computation(geometry)) + generator = HacklGenerator_OneLambda(dummy, R_STATOR_END, offset=OFFSET) + + np.random.seed(SEED) + while True: + generator.set_parameters(generator.random_parameters()) + barriers = generator.split_barriers(generator.generate_barriers()) + if generator.feasible_barriers(barriers): + return barriers + + +@pytest.fixture(scope="module") +def solved(tmp_path_factory): + """One solve, shared by both subjects -- an FEA solve is far too costly to repeat.""" + geometry = Geometry() + computation = ComputationWithPhaseQuantities(geometry) + barriers = _feasible_barriers(geometry) + + design = None + try: + design = Design.create( + "TransformBoundary", + "Design01", + str(tmp_path_factory.mktemp("aedt") / "transform.aedt"), + geometry, + computation, + version=AEDT_VERSION, + non_graphical=True, + new_desktop=False, + close_on_exit=False, + ) + design.add_rotor() + for barrier in barriers: + design.add_rotor_barrier(barrier) + return design.compute(*CURRENTS, NUM_CORES=NUM_CORES) + finally: + if design is not None: + design.close_project() + + +@pytest.mark.parametrize(("label", "phase_exprs", "prefix"), SUBJECTS, ids=[s[0] for s in SUBJECTS]) +@pytest.mark.parametrize("harmonic", [1, 3]) +def test_python_transform_matches_ansys(solved, label, phase_exprs, prefix, harmonic): + assert solved is not None, "solve returned no results" + needed = [*phase_exprs, "theta_el", f"{prefix}_d{harmonic}", f"{prefix}_q{harmonic}"] + missing = [k for k in needed if k not in solved] + assert not missing, f"missing from the extracted solution: {missing}" + + # (n_timesteps, 5), phases in order A..E; to_dq reads the phase count off that axis. + per_phase = np.column_stack([np.asarray(solved[k], dtype=float) for k in phase_exprs]) + # Ansys returns theta_el in DEGREES (its "- pi" evaluates as -180 in an angular + # context); transforms.to_dq takes radians. + theta_el = np.deg2rad(np.asarray(solved["theta_el"], dtype=float)) + + d_py, q_py = to_dq(per_phase, theta_el, harmonic) + d_aedt = np.asarray(solved[f"{prefix}_d{harmonic}"], dtype=float) + q_aedt = np.asarray(solved[f"{prefix}_q{harmonic}"], dtype=float) + + # rtol handles solver/report round-off; atol keeps near-zero samples from failing on + # relative error alone, scaled to the signal rather than to absolute units. + atol = 1e-6 * max(np.max(np.abs(d_aedt)), np.max(np.abs(q_aedt)), 1.0) + np.testing.assert_allclose(d_py, d_aedt, rtol=1e-6, atol=atol, err_msg=f"{label} d{harmonic} mismatch") + np.testing.assert_allclose(q_py, q_aedt, rtol=1e-6, atol=atol, err_msg=f"{label} q{harmonic} mismatch") diff --git a/tests/test_synrm_5f_40s.py b/tests/test_synrm_5f_40s.py index 9d077b1..75fa0b9 100644 --- a/tests/test_synrm_5f_40s.py +++ b/tests/test_synrm_5f_40s.py @@ -3,6 +3,7 @@ import numpy as np import pytest +from machine_design.transforms import electrical_angle, to_dq from motors.synrm_3f_36s import Geometry as BaseGeometry from motors.synrm_5f_40s import Computation, Geometry @@ -67,8 +68,10 @@ def test_derived_params_are_set(geometry2): def test_solution_expressions_is_a_list_ending_in_torque(computation2): assert isinstance(computation2.solution_expressions, list) - assert computation2.solution_expressions[-1] == "Moving1.Torque" - assert len(computation2.solution_expressions) == 23 + assert computation2.solution_expressions[0] == "Moving1.Position" + assert computation2.solution_expressions[1] == "Moving1.Torque" + assert computation2.solution_expressions[-1] == "L(PhaseE,PhaseE)" + assert len(computation2.solution_expressions) == 42 def test_post_params_plot_names_unique(computation2): @@ -76,14 +79,29 @@ def test_post_params_plot_names_unique(computation2): assert len(plot_names) == len(set(plot_names)) -def test_output_vars_known_values(computation2): - assert computation2.output_vars["Rstat"] == "19" - assert computation2.output_vars["Lew"] == "0" +def test_output_vars_is_empty(computation2): + assert computation2.output_vars == {} -def test_extract_results_converts_units_by_prefix(computation2): - solutions = SimpleNamespace(data_real=lambda expr: 100.0) +def test_Rstat_Lew_values(computation2): + assert computation2.Rstat == 19.0 + assert computation2.Lew == 0.0 + + +def test_extract_results_converts_units(computation2): + computation2.set_variables(SimpleNamespace(variable_manager={}), Id1=1.0, Iq1=2.0, Id3=3.0, Iq3=4.0) + + def data_real(expr): + letters = ["ABCDE".index(part[0]) + 1 for part in expr.split("Phase")[1:]] + value = float(np.prod(letters) if letters else 7.0) + return np.full(2, value * 1e9 if expr.startswith("L(") else value) + + solutions = SimpleNamespace(data_real=data_real, primary_sweep_values=[0.0, 1.0]) out = computation2.extract_results(solutions) - assert out["Ld1"] == pytest.approx(100.0 / 1e9) # Ld*/Lq* , nH - assert out["I_d1"] == pytest.approx(100.0 / 1e3) # I_* , A - assert out["V_d1"] == pytest.approx(100.0) + + theta_el = electrical_angle(np.full(2, 7.0), computation2.InitPos, computation2.geometry.PolePairs, computation2.RotSign, degrees=True) + current_phases = np.stack([np.full(2, i + 1.0) for i in range(5)], axis=-1) + expected_I_d1, _ = to_dq(current_phases, theta_el, harmonic=1) + + assert out["I_d1"] == pytest.approx(expected_I_d1 / 1e3) + assert out["Moving1.Torque"] == pytest.approx(np.full(2, 7.0)) diff --git a/tests/test_transforms.py b/tests/test_transforms.py new file mode 100644 index 0000000..6fc8fd8 --- /dev/null +++ b/tests/test_transforms.py @@ -0,0 +1,128 @@ +"""Pins the m-phase dq convention used by `machine_design.transforms`. + +These need no Ansys. They exist because the same transform is written out by hand in +`motors/synrm_5f_40s.py`'s output variables, and a silent disagreement between the two +would corrupt every voltage and flux number downstream without failing loudly. The live +cross-check against Ansys is `tests/test_ansys_transforms.py`. + +Parametrised over phase count as well as harmonic, so the 3-phase anchor is covered by +the same transform as the 5-phase one. +""" + +import numpy as np +import pytest + +from machine_design.transforms import electrical_angle, to_dq, to_phases + +THETA = np.linspace(0.0, 2.0 * np.pi, 37) + +# (n_phases, harmonic). A harmonic that is a multiple of the phase count is +# zero-sequence and not invertible -- see test_triplen_harmonic_is_zero_sequence. +CASES = [(3, 1), (5, 1), (5, 3)] +CASE_IDS = [f"{m}ph-h{h}" for m, h in CASES] + + +@pytest.mark.parametrize(("n_phases", "harmonic"), CASES, ids=CASE_IDS) +@pytest.mark.parametrize(("d", "q"), [(1.0, 0.0), (0.0, 1.0), (0.7, -0.3), (-1.2, 2.5)]) +def test_round_trip_recovers_dq(n_phases, harmonic, d, q): + x = to_phases(d, q, THETA, harmonic, n_phases) + assert x.shape == (len(THETA), n_phases) + d_out, q_out = to_dq(x, THETA, harmonic) + np.testing.assert_allclose(d_out, d, atol=1e-12) + np.testing.assert_allclose(q_out, q, atol=1e-12) + + +@pytest.mark.parametrize(("h_signal", "h_read"), [(1, 3), (3, 1)]) +def test_harmonics_are_orthogonal(h_signal, h_read): + """A pure dq1 waveform must read as zero in dq3, and vice versa. + + This is what makes the two current subspaces independent design variables; if it + broke, dq3 injection would leak into the dq1 measurement. + """ + x = to_phases(1.3, -0.8, THETA, h_signal, 5) + d, q = to_dq(x, THETA, h_read) + np.testing.assert_allclose(d, 0.0, atol=1e-12) + np.testing.assert_allclose(q, 0.0, atol=1e-12) + + +def test_superposed_harmonics_separate(): + """dq1 and dq3 components ride on the same phase waveform without interfering.""" + x = to_phases(2.0, 0.5, THETA, 1, 5) + to_phases(-0.4, 0.9, THETA, 3, 5) + d1, q1 = to_dq(x, THETA, 1) + d3, q3 = to_dq(x, THETA, 3) + np.testing.assert_allclose(d1, 2.0, atol=1e-12) + np.testing.assert_allclose(q1, 0.5, atol=1e-12) + np.testing.assert_allclose(d3, -0.4, atol=1e-12) + np.testing.assert_allclose(q3, 0.9, atol=1e-12) + + +def test_triplen_harmonic_is_zero_sequence_in_three_phase(): + """m = 3, h = 3 is deliberately NOT invertible, and that is physics. + + 3*(2*pi*k/3) is a whole number of turns, so every phase sees the same angle: the + triplen harmonic is zero-sequence. Pinned so nobody "fixes" the transform to make + this round-trip. + """ + x = to_phases(1.0, 0.0, THETA, 3, 3) + np.testing.assert_allclose(x[:, 0], x[:, 1], atol=1e-12) + np.testing.assert_allclose(x[:, 0], x[:, 2], atol=1e-12) + + d, _ = to_dq(x, THETA, 3) + assert not np.allclose(d, 1.0, atol=1e-6), "expected the triplen round-trip to fail" + + +@pytest.mark.parametrize("n_phases", [3, 5]) +def test_amplitude_scaling_is_two_over_m(n_phases): + """The 2/m factor: a unit-amplitude d-axis waveform must give d = 1.""" + d, _ = to_dq(to_phases(1.0, 0.0, THETA, 1, n_phases), THETA, 1) + np.testing.assert_allclose(d, 1.0, atol=1e-12) + + +@pytest.mark.parametrize("n_phases", [3, 5]) +def test_q_axis_sign_convention(n_phases): + """`q` is the NEGATIVE sine sum, matching sin(-h*(theta - 2*pi*k/m)) in Ansys. + + Built from explicit cos/sin phase quantities rather than via `to_phases`, so this + fails if the sign in `to_dq` is ever "corrected" to the more usual convention. + """ + theta = 0.37 + phi = theta - 2.0 * np.pi * np.arange(n_phases) / n_phases + + # A pure +cos(phi) phase quantity is pure d-axis, unit amplitude. + d, q = to_dq(np.cos(phi), theta, 1) + assert d == pytest.approx(1.0, abs=1e-12) + assert q == pytest.approx(0.0, abs=1e-12) + + # A pure +sin(phi) phase quantity gives q = -1 under this convention, not +1. + d, q = to_dq(np.sin(phi), theta, 1) + assert d == pytest.approx(0.0, abs=1e-12) + assert q == pytest.approx(-1.0, abs=1e-12) + + +def test_electrical_angle_matches_ansys_expression(): + """theta_el = RotSign*(Position - InitPos)*PolePairs - pi, in radians.""" + position = np.array([0.0, 0.25, 0.5]) + init_pos = np.deg2rad(-45.0) + got = electrical_angle(position, init_pos, pole_pairs=2, rot_sign=1.0) + expected = 1.0 * (position - init_pos) * 2 - np.pi + np.testing.assert_allclose(got, expected, atol=1e-12) + + +def test_electrical_angle_accepts_degrees_like_ansys_reports_them(): + """Ansys reports `Moving1.Position` in degrees; `degrees=True` must match radians. + + Guards the trap this API invites: the natural input comes back from Ansys in + degrees, and passing it as radians fails silently rather than loudly. + """ + position_deg, init_deg = np.array([315.0, 316.8, 318.6]), -45.0 + from_deg = electrical_angle(position_deg, init_deg, pole_pairs=2, degrees=True) + from_rad = electrical_angle(np.deg2rad(position_deg), np.deg2rad(init_deg), pole_pairs=2) + np.testing.assert_allclose(from_deg, from_rad, atol=1e-12) + + # The observed Ansys value: (315 - -45) * 2 - 180 degrees = 540 deg = 3*pi rad. + assert from_deg[0] == pytest.approx(3.0 * np.pi, abs=1e-12) + + +def test_rejects_too_few_phases(): + with pytest.raises(ValueError, match="at least 3 phases"): + to_dq(np.zeros((len(THETA), 2)), THETA, 1)