Bugfix in BEM hydro forces: using the right incidences now

When transforming excitation coefficients to be relative to the incident wave heading rather than to the global frame, we were using the unsorted headings vector (heads) instead of the sorted headings vector (self.BEM_headings)

raft/raft_fowt.py

@@ -665,13 +665,15 @@ def calcBEM(self, dw=0, wMax=0, wInf=10.0, dz=0, da=0, headings=[0], meshDir=os.
         # The Tflag means that the first column is in units of periods, not frequencies, and therefore the first set (-1) becomes zero-frequency and the second set is infinite
 
         # process and sort headings and sort frequencies
+        R_unsorted = R.copy()
+        I_unsorted = I.copy()
         self.BEM_headings = np.array(heads)%(360)  # save headings in range of 0-360 [deg]            # interpole to the frequencies RAFT is using        
         sorted_indices = np.argsort(self.BEM_headings)
         self.BEM_headings = self.BEM_headings[sorted_indices]
         M = M[sorted_indices,:,:]
         P = P[sorted_indices,:,:]
         R = R[sorted_indices,:,:]
-        I = I[sorted_indices,:,:]            
+        I = I[sorted_indices,:,:]
 
         # interpolate to RAFT model frequencies
         # zero frequency values are being stacked on to give smooth results if the requested frequency is below what's available from HAMS
@@ -692,9 +694,9 @@ def calcBEM(self, dw=0, wMax=0, wInf=10.0, dz=0, da=0, headings=[0], meshDir=os.
         # for accurate magnitudes when interpolating between directions.
         self.X_BEM = np.zeros_like(X_BEM_temp)
 
-        for ih in range(len(heads)):
-            sin_heading = np.sin(np.radians(heads[ih]))
-            cos_heading = np.cos(np.radians(heads[ih]))
+        for ih in range(len(self.BEM_headings)):
+            sin_heading = np.sin(np.radians(self.BEM_headings[ih]))
+            cos_heading = np.cos(np.radians(self.BEM_headings[ih]))
 
             self.X_BEM[ih,0,:] =  cos_heading * X_BEM_temp[ih,0,:] + sin_heading * X_BEM_temp[ih,1,:]
             self.X_BEM[ih,1,:] = -sin_heading * X_BEM_temp[ih,0,:] + cos_heading * X_BEM_temp[ih,1,:]