Adding 2D plotting ability for rotors:

- Changes in a few levels of plotting commands to support
  2D plotting of rotors including options to show the blades
  or just draw a line/circle for a more simple view.
- A few usage things probably need cleaning up to support all
  options. I left drawing a tip circle on by accident...

raft/raft_fowt.py

@@ -2081,11 +2081,21 @@ def plot(self, ax, color=None, nodes=0, plot_rotor=True, station_plot=[],
         # including hydro excitation vectors stored in each member
 
 
-    def plot2d(self, ax, color=None, station_plot=[], Xuvec=[1,0,0], Yuvec=[0,0,1]):
-        '''plots the FOWT in 2d - only the platform and moorings so far...'''
+    def plot2d(self, ax, color=None, plot_rotor=1, 
+        Xuvec=[1,0,0], Yuvec=[0,0,1]):
+        '''plots the FOWT in 22.
+        
+        Parameters
+        ----------
+        plot_rotor : int
+            0 : don't plot; 1: plot 3D view; 2: plot a simple line
+        '''
 
         R = rotationMatrix(self.r6[3], self.r6[4], self.r6[5])  # note: eventually Rotor could handle orientation internally <<<
-
+
+        Xuvec = np.array(Xuvec)
+        Yuvec = np.array(Yuvec)
+
         if self.ms:
             self.ms.plot2d(ax=ax, color=color, Xuvec=Xuvec, Yuvec=Yuvec)
 
@@ -2096,3 +2106,40 @@ def plot2d(self, ax, color=None, station_plot=[], Xuvec=[1,0,0], Yuvec=[0,0,1]):
         for mem in self.memberList:
             mem.setPosition()
             mem.plot(ax, r_ptfm=self.r6[:3], R_ptfm=R, color=color, plot2d=True, Xuvec=Xuvec, Yuvec=Yuvec)
+
+        # Plot the rotor(s)
+        if plot_rotor == 1:
+            for rotor in self.rotorList:
+                coords = np.array([rotor.coords[0], rotor.coords[1], 0]) + self.r6[:3]
+                rotor.plot(ax, r_ptfm=coords, R_ptfm=R, color=color, plot2d=True,
+                Xuvec=Xuvec, Yuvec=Yuvec)
+
+        elif plot_rotor == 2:
+            # simple circle/line plot of rotor, ignoring precone or tilt
+            for rotor in self.rotorList:
+                r = rotor.ccblade.r[-1]  # rotor tip radius [m]
+
+                X=[]  # lists to be filled with coordinates for plotting
+                Y=[]
+                Z=[]
+
+                n = 24  # number of sides for a circle
+                for i in range(n+1):
+                    y = np.cos(float(i)/float(n)*2.0*np.pi)    # x coordinates of a unit circle
+                    z = np.sin(float(i)/float(n)*2.0*np.pi)    # y
+
+                    X.append(0)
+                    Y.append(r*y)
+                    Z.append(r*z)
+
+                R_heading = rotationMatrix(0, 0, rotor.heading)  # rotation matrix for rotor heading
+
+                P2 = np.vstack([X, Y, Z])  # combine into a matrix of coordinates
+                #P2 = P2 + np.array([-rotor.overhang, 0, rotor.Zhub])[:,None]  # overhang and height
+                #P2 = np.matmul(np.matmul(R, R_heading), P2) + rotor.r3[:,None]  # with ptm rotation
+                P2 = np.matmul(R_heading, P2) + rotor.r3[:,None]  # with only heading
+
+                # apply any 3D to 2D transformation here to provide desired viewing angle
+                Xs2d = np.matmul(Xuvec, P2)
+                Ys2d = np.matmul(Yuvec, P2)
+                ax.plot(Xs2d, Ys2d, color=color, lw=1.0)

raft/raft_model.py

@@ -1510,7 +1510,8 @@ def plot(self, ax=None, hideGrid=False, draw_body=True, color=None, nodes=0,
 
 
     def plot2d(self, ax=None, hideGrid=False, draw_body=True, color=None, 
-               station_plot=[], Xuvec=[1,0,0], Yuvec=[0,0,1], figsize=(6,4)):
+               station_plot=[], Xuvec=[1,0,0], Yuvec=[0,0,1], figsize=(6,4),
+               plot_rotor=2):
         '''plots the whole model, including FOWTs and mooring system...'''
 
         # for now, start the plot via the mooring system, since MoorPy doesn't yet know how to draw on other codes' plots
@@ -1534,7 +1535,7 @@ def plot2d(self, ax=None, hideGrid=False, draw_body=True, color=None,
 
         # plot each FOWT
         for fowt in self.fowtList:
-            fowt.plot2d(ax, color=color, station_plot=station_plot, Xuvec=Xuvec, Yuvec=Yuvec)
+            fowt.plot2d(ax, color=color, plot_rotor=plot_rotor, Xuvec=Xuvec, Yuvec=Yuvec)
 
         ax.axis("equal")
 

raft/raft_rotor.py

@@ -50,6 +50,9 @@ def __init__(self, turbine, w, ir):
 
         self.nBlades    = getFromDict(turbine, 'nBlades', shape=turbine['nrotors'], dtype=int)[ir]         # [-]
 
+        self.heading = 0  # heading of the rotor in radians
+
+        # note: the below are not headings, they're blade azimuth angles -- need to rename
         default_headings    = list(np.arange(self.nBlades) * 360 / self.nBlades) # equally distribute blades
         self.headings       = getFromDict(turbine, 'headings', shape=-1, default=default_headings)  # [deg]
 
@@ -782,6 +785,8 @@ def calcAero(self, case, current=False, display=0):
         # inflow misalignment heading relative to turbine heading [deg]
         yaw_misalign =  turbine_heading - np.radians(heading)
         turbine_tilt = turbine_tilt
+
+        self.heading = turbine_heading  # store for later use (since q doesn't yet account for nacelle yaw!)
 
         # call CCBlade
         loads, derivs = self.runCCBlade(speed, ptfm_pitch=turbine_tilt,
@@ -1000,9 +1005,20 @@ def calcAero(self, case, current=False, display=0):
 
 
     def plot(self, ax, r_ptfm=[0,0,0], R_ptfm=np.eye(3), azimuth=0, color='k', 
-             airfoils=False, zorder=2):
-        '''Draws the rotor on the passed axes, considering optional platform offset and rotation matrix, and rotor azimuth angle'''
+             airfoils=False, plot2d=False, Xuvec=[1,0,0], Yuvec=[0,0,1], zorder=2):
+        '''Draws the rotor on the passed axes, considering optional platform 
+        offset and rotation matrix, and rotor azimuth angle.
+        
+        Parameters
+        ----------        
+        plot2d: bool
+            If true, produces a 2d plot on the axes defined by Xuvec and Yuvec. 
+            Otherwise produces a 3d plot (default).
+        '''
 
+        Xuvec = np.array(Xuvec)
+        Yuvec = np.array(Yuvec)
+
         # support self color option
         if color == 'self':
             color = self.color  # attempt to allow custom colors
@@ -1034,7 +1050,7 @@ def plot(self, ax, r_ptfm=[0,0,0], R_ptfm=np.eye(3), azimuth=0, color='k',
 
         P = np.array([X, Y, Z])
 
-        # ----- rotation matricse ----- 
+        # ----- rotation matrices ----- 
         # (blade pitch would be a -rotation about local z)
         R_precone = rotationMatrix(0, -self.ccblade.precone, 0)  
         #R_azimuth = [rotationMatrix(azimuth + azi, 0, 0) for azi in 2*np.pi/3.*np.arange(3)]
@@ -1050,20 +1066,64 @@ def plot(self, ax, r_ptfm=[0,0,0], R_ptfm=np.eye(3), azimuth=0, color='k',
             P2 = P2 + np.array([-self.overhang, 0, self.Zhub])[:,None] # PRP to tower-shaft intersection point      # assumes overhang value is always positive
             P2 = np.matmul(R_ptfm, P2) + np.array(r_ptfm)[:,None]
 
-            # drawing airfoils                            
-            if airfoils:
-                for ii in range(m-1):
-                    ax.plot(P2[0, npts*ii:npts*(ii+1)], P2[1, npts*ii:npts*(ii+1)], P2[2, npts*ii:npts*(ii+1)], color=color, lw=0.4)  
-            # draw outline
-            ax.plot(P2[0, 0:-1:npts], P2[1, 0:-1:npts], P2[2, 0:-1:npts], color=color, lw=0.4, zorder=zorder) # leading edge  
-            ax.plot(P2[0, 2:-1:npts], P2[1, 2:-1:npts], P2[2, 2:-1:npts], color=color, lw=0.4, zorder=zorder)  # trailing edge
-
+            if plot2d:  # new 2d plotting option
+
+                # apply any 3D to 2D transformation here to provide desired viewing angle
+                Xs2d = np.matmul(Xuvec, P2)
+                Ys2d = np.matmul(Yuvec, P2)
+
+                if airfoils:
+                    for ii in range(m-1):
+                        ax.plot(Xs2d[npts*ii:npts*(ii+1)], Ys2d[npts*ii:npts*(ii+1)], color=color, lw=0.4)  
+                # draw outline
+                ax.plot(Xs2d[0:-1:npts], Ys2d[0:-1:npts], color=color, lw=0.4, zorder=zorder) # leading edge  
+                ax.plot(Xs2d[2:-1:npts], Ys2d[2:-1:npts], color=color, lw=0.4, zorder=zorder)  # trailing edge
+
 
-        #for j in range(m):
-        #    linebit.append(ax.plot(Xs[j::m], Ys[j::m], Zs[j::m]            , color='k'))  # station rings
-        #
-        #return linebit
+            else:  # normal 3d case
+                # drawing airfoils                            
+                if airfoils:
+                    for ii in range(m-1):
+                        ax.plot(P2[0, npts*ii:npts*(ii+1)], P2[1, npts*ii:npts*(ii+1)], P2[2, npts*ii:npts*(ii+1)], color=color, lw=0.4)  
+                # draw outline
+                ax.plot(P2[0, 0:-1:npts], P2[1, 0:-1:npts], P2[2, 0:-1:npts], color=color, lw=0.4, zorder=zorder) # leading edge  
+                ax.plot(P2[0, 2:-1:npts], P2[1, 2:-1:npts], P2[2, 2:-1:npts], color=color, lw=0.4, zorder=zorder)  # trailing edge
+
+        # ----- Also draw a circle -----
+        r = self.ccblade.r[-1]
+        x_shift = r * np.tan(self.ccblade.precone)  # shift forward of tips due to precone
+
+        # lists to be filled with coordinates for plotting
+        X=[]
+        Y=[]
+        Z=[]
+
+        n = 24                      # number of sides for a circle
+        for i in range(n+1):
+            y = np.cos(float(i)/float(n)*2.0*np.pi)    # x coordinates of a unit circle
+            z = np.sin(float(i)/float(n)*2.0*np.pi)    # y
+
+            X.append(0)
+            Y.append(r*y)
+            Z.append(r*z)
 
+        Pcirc = np.vstack([X, Y, Z])
+        P2 = np.matmul(R_tilt, Pcirc)
+        # PRP to tower-shaft intersection point 
+        # Assumes overhang value is always positive, also accounts for precurve
+        P2 = P2 + np.array([-self.overhang - x_shift, 0, self.Zhub])[:,None] 
+        P2 = np.matmul(R_ptfm, P2) + np.array(r_ptfm)[:,None]
+
+        if plot2d:  # new 2d plotting option
+            # apply any 3D to 2D transformation here to provide desired viewing angle
+            Xs2d = np.matmul(Xuvec, P2)
+            Ys2d = np.matmul(Yuvec, P2)
+            ax.plot(Xs2d, Ys2d, color=color, lw=0.4, zorder=zorder)
+
+        else:  # normal 3d case
+            ax.plot(P2[0,:], P2[1,:], P2[2,:], color=color, lw=0.4, zorder=zorder)
+
+
 
     def IECKaimal(self, case, current=False):        # 
         '''Calculates rotor-averaged turbulent wind spectrum based on inputted turbulence intensity or class.'''