From efa515f2895c62d204c260ea02e76c62858f41c0 Mon Sep 17 00:00:00 2001
From: James Packer <james.packer@ansys.com>
Date: Fri, 13 Dec 2024 10:18:56 +0000
Subject: [PATCH 1/2] Initial implementation for graph harmonics

---
 .../core/methods/rpc_methods_graphs.py        | 63 +++++++++++++++++++
 1 file changed, 63 insertions(+)

diff --git a/src/ansys/motorcad/core/methods/rpc_methods_graphs.py b/src/ansys/motorcad/core/methods/rpc_methods_graphs.py
index 4931d3775..7ac5feb16 100644
--- a/src/ansys/motorcad/core/methods/rpc_methods_graphs.py
+++ b/src/ansys/motorcad/core/methods/rpc_methods_graphs.py
@@ -22,6 +22,9 @@
 
 """RPC methods for graphs."""
 from dataclasses import dataclass
+import math
+
+import numpy
 
 from ansys.motorcad.core.rpc_client_core import MotorCADError
 
@@ -220,6 +223,66 @@ def get_magnetic_graph(self, graph_name):
         else:
             return self._get_graph(self.get_magnetic_graph_point, graph_name)
 
+    def get_magnetic_graph_harmonics(self, graph_name):
+        """Get harmonic analysis from Motor-CAD magnetic graph.
+
+        Parameters
+        ----------
+        graph_name : str, int
+            Name (preferred) or ID of the graph. In Motor-CAD, you can
+            select **Help -> Graph Viewer** to see the graph name.
+        Returns
+        -------
+        order_values : list
+            Value of harmonic orders from graph
+        amplitude_values : list
+            Value of harmonic amplitudes from graph
+        angle_values : list
+            Value of harmonic angles from graph in degrees
+        """
+        x, y = self.get_magnetic_graph(graph_name)
+        # Find x-axis limits and range, as this is needed to find the phase information
+        min_x = min(x)
+        cycles = (max(x) - min(x)) / 360
+
+        # y normally contains a duplicated final point (360deg=0deg), so discard this point.
+        y_no_duplicate = y[: len(y) - 1]
+
+        # Carry out FFT only get up to the Nyquist limit, using real valued inputs
+        y_fft = numpy.fft.rfft(y_no_duplicate, norm="forward")
+
+        # Apply normalisation.
+        for i in range(len(y_fft)):
+            if i > 0:
+                # Multiply by 2 to account for the positive and negative frequency component
+                y_fft[i] = 2 * y_fft[i]
+            else:
+                # 0th component does not need to be doubled
+                y_fft[i] = y_fft[i]
+
+        # Get amplitude and angle:
+        y_mag = numpy.absolute(y_fft)
+        y_ang = numpy.angle(y_fft, deg=True)
+
+        # Motor-CAD harmonic plot convention shifts the angles by 90 degrees.
+        # Also consider the phase angle of the first point:
+        for i in range(len(y_ang)):
+            y_ang[i] = y_ang[i] + 90 - (min_x * i / cycles)
+            if y_ang[i] > 180:
+                y_ang[i] = y_ang[i] - 360
+            if y_ang[i] < -180:
+                y_ang[i] = y_ang[i] + 360
+
+        # For very small magnitudes, the angle information is not meaningful, so set to zero
+        for i in range(len(y_ang)):
+            if math.isclose(y_mag[i], 0, abs_tol=1e-8):
+                y_ang[i] = 0
+
+        # Generate an index list for plotting
+        y_index = numpy.linspace(start=0, stop=len(y_fft) / cycles, num=len(y_fft), endpoint=False)
+
+        return [y_index.tolist(), y_mag.tolist(), y_ang.tolist()]
+
     def get_temperature_graph(self, graph_name):
         """Get graph points from a Motor-CAD transient temperature graph.
 

From b3ecc74ce0b1ec48b4d30efa23f8a585dff7ab9b Mon Sep 17 00:00:00 2001
From: James Packer <james.packer@ansys.com>
Date: Fri, 13 Dec 2024 10:37:39 +0000
Subject: [PATCH 2/2] Add test

---
 tests/test_graphs.py | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/tests/test_graphs.py b/tests/test_graphs.py
index f66511532..e6e343c9c 100644
--- a/tests/test_graphs.py
+++ b/tests/test_graphs.py
@@ -169,6 +169,24 @@ def test_get_magnetic_3d_graph(mc):
     assert almost_equal(graph_result.data[1][1], y)
 
 
+def test_get_magnetic_graph_harmonics(mc):
+    reset_to_default_file(mc)
+
+    mc.set_variable("TorqueCalculation", True)
+    mc.set_variable("ElectromagneticForcesCalc_Load", True)
+
+    mc.do_magnetic_calculation()
+    orders, amp, ang = mc.get_magnetic_graph_harmonics("PhaseCurrent1")
+    amplitude_expected = mc.get_variable("PeakCurrent")
+    phase_expected = mc.get_variable("PhaseAdvance") - 90
+    assert len(orders) == 2
+    assert len(amp) == 2
+    assert len(ang) == 2
+    assert orders[1] == 1
+    assert almost_equal(amp[1], amplitude_expected)
+    assert almost_equal(ang[1], phase_expected)
+
+
 #   #Not fully ready submitted an issue
 # def test_get_fea_graph_point():
 #     reset_to_default_file(mc)