Merge branch 'new_pc_file' into 'master'

New pc file

See merge request toolbox/WindEnergyToolbox!43

wetb/hawc2/bladeinfo.py → wetb/hawc2/blade.py

@@ -11,6 +11,8 @@ from scipy.interpolate.interpolate import interp1d
 import os
 from wetb.utils.geometry import rad
 from wetb.hawc2.st_file import StFile
+from wetb.hawc2.ae_file import AEFile
+from wetb.hawc2.mainbody import MainBody
 
 
 
@@ -61,7 +63,7 @@ from wetb.hawc2.st_file import StFile
 #         """
 
 
-class H2aeroBladeInfo(PCFile, AtTimeFile):
+class H2aeroBlade(PCFile, AEFile, AtTimeFile):
     """Provide HAWC2 info about a blade
     
     From AE file:
@@ -86,36 +88,37 @@ class H2aeroBladeInfo(PCFile, AtTimeFile):
         
     """
 
-    def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, blade_name=None):
-        
-        if isinstance(htcfile, str):
-            assert htcfile.lower().endswith('.htc')
-            htcfile = HTCFile(htcfile)
-        self.htcfile = htcfile
-        blade_name = blade_name or htcfile.aero.link[2]
-        #s = htcfile.new_htc_structure
-        at_time_filename = at_time_filename or ("output_at_time" in htcfile and os.path.join(htcfile.modelpath, htcfile.output_at_time.filename[0] + ".dat"))
-        pc_filename = pc_filename or os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0])
-        ae_filename = ae_filename or os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0])
-        
+    def __init__(self, htcfile=None, ae_filename=None, pc_filename=None, at_time_filename=None, blade_name=None):
+        if htcfile:
+            if isinstance(htcfile, str):
+                assert htcfile.lower().endswith('.htc')
+                htcfile = HTCFile(htcfile)
+            self.htcfile = htcfile
+            blade_name = blade_name or htcfile.aero.link[2]
+            at_time_filename = at_time_filename or ("output_at_time" in htcfile and os.path.join(htcfile.modelpath, htcfile.output_at_time.filename[0] + ".dat"))
+            pc_filename = pc_filename or os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0])
+            ae_filename = ae_filename or os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0])
+            self.hawc2_splines_data = self.hawc2_splines()
+            
         #mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
         #self.mainbody_blade = [mb for mb in mainbodies if mb.name[0] == blade_name][0]
         
         if os.path.isfile(pc_filename) and os.path.isfile(ae_filename):
-            PCFile.__init__(self, pc_filename, ae_filename)
+            AEFile.__init__(self, ae_filename)
+            PCFile.__init__(self, pc_filename)
             blade_radius = self.ae_sets[1][-1,0]
         
-        if os.path.isfile(at_time_filename):
+        if at_time_filename and os.path.isfile(at_time_filename):
             AtTimeFile.__init__(self, at_time_filename, blade_radius)
             self.curved_length = self.radius_curved_ac()[-1]
         else:
-            raise NotImplementedError
+            self.curved_length = None
             #z_nd = (np.cos(np.linspace(np.pi, np.pi*2,len(curved_length)-1))+1)/2
             #self.curved_length = np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))[-1]
 
         
         
-        self.hawc2_splines_data = self.hawc2_splines()
+        
         
         
     @property
@@ -222,8 +225,103 @@ class H2aeroBladeInfo(PCFile, AtTimeFile):
                 assert np.all(l>=self.c2def[0,2]) and np.all(l<=self.c2def[-1,2])
                 return self.c2nd(l/self.c2def[-1, 2])
             
+    def _Cxxx(self, radius, alpha, column, ae_set_nr=1):
+        thickness = self.thickness(radius, ae_set_nr)
+        pc_set_nr = self.pc_set_nr(radius, ae_set_nr)
+        thicknesses, profiles = self.pc_sets[pc_set_nr]
+        index = np.searchsorted(thicknesses, thickness)
+        if index == 0:
+            index = 1
+
+        Cx0, Cx1 = profiles[index - 1:index + 1]
+        Cx0 = np.interp(alpha, Cx0[:, 0], Cx0[:, column])
+        Cx1 = np.interp(alpha, Cx1[:, 0], Cx1[:, column])
+        th0, th1 = thicknesses[index - 1:index + 1]
+        return Cx0 + (Cx1 - Cx0) * (thickness - th0) / (th1 - th0)
+    
+    def _CxxxH2(self, radius, alpha, column, ae_set_nr=1):
+        thickness = self.thickness(radius, ae_set_nr)
+        pc_set_nr = self.pc_set_nr(radius, ae_set_nr)
+        thicknesses, profiles = self.pc_sets[pc_set_nr]
+        index = np.searchsorted(thicknesses, thickness)
+        if index == 0:
+            index = 1
+
+        Cx0, Cx1 = profiles[index - 1:index + 1]
+        
+        Cx0 = np.interp(np.arange(360), Cx0[:,0]+180, Cx0[:,column])
+        Cx1 = np.interp(np.arange(360), Cx1[:,0]+180, Cx1[:,column])
+        #Cx0 = np.interp(alpha, Cx0[:, 0], Cx0[:, column])
+        #Cx1 = np.interp(alpha, Cx1[:, 0], Cx1[:, column])
+        th0, th1 = thicknesses[index - 1:index + 1]
+        cx = Cx0 + (Cx1 - Cx0) * (thickness - th0) / (th1 - th0)
+        return np.interp(alpha+180, np.arange(360), cx)
+    
+        
+
+    def CL(self, radius, alpha, ae_set_nr=1):
+        """Lift coefficient
+
+        Parameters
+        ---------
+        radius : float
+            radius [m]
+        alpha : float
+            Angle of attack [deg]
+        ae_set_nr : int optional
+            Aerdynamic set number, default is 1
+
+        Returns
+        -------
+        Lift coefficient : float
+        """
+        return self._Cxxx(radius, alpha, 1, ae_set_nr)
+
+
+    def CL_H2(self, radius, alpha, ae_set_nr=1):
+        return self._CxxxH2(radius, alpha, 1, ae_set_nr)
+    
+    def CD(self, radius, alpha, ae_set_nr=1):
+        """Drag coefficient
+
+        Parameters
+        ---------
+        radius : float
+            radius [m]
+        alpha : float
+            Angle of attack [deg]
+        ae_set_nr : int optional
+            Aerdynamic set number, default is 1
+
+        Returns
+        -------
+        Drag coefficient : float
+        """
+        return self._Cxxx(radius, alpha, 2, ae_set_nr)
+
+    def CM(self, radius, alpha, ae_set_nr=1):
+        return self._Cxxx(radius, alpha, 3, ae_set_nr)
             
-class H2BladeInfo(H2aeroBladeInfo):
+            
+    def plot_xz_geometry(self, plt):
+
+        z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]), label='Center line')
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) + self.chord(z) / 2, label='Leading edge')
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) - self.chord(z) / 2, label="Trailing edge")
+        curve_l_nd, x, y, z, twist = self.hawc2_splines()
+        plt.plot(z, x, label='Hawc2spline')
+
+    def plot_yz_geometry(self, plt):
+
+        z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]), label='Center line')
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.thickness(z) / 100 * self.chord(z) / 2, label='Suction side')
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) - self.thickness(z) / 100 * self.chord(z) / 2, label="Pressure side")
+        curve_l_nd, x, y, z, twist = self.hawc2_splines()
+        plt.plot(z, y, label='Hawc2spline')
+            
+class H2Blade(H2aeroBlade, MainBody):
     """Provide HAWC2 info about a blade
     
     From AE file:
@@ -247,19 +345,24 @@ class H2BladeInfo(H2aeroBladeInfo):
     with template
         
     """
-    def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
-        if isinstance(htcfile, str):
-            htcfile = HTCFile(htcfile)
-        s = htcfile.new_htc_structure
+    def __init__(self, htcfile=None, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
+        if htcfile is not None:
+            if isinstance(htcfile, str):
+                htcfile = HTCFile(htcfile)
+            s = htcfile.new_htc_structure
 #         at_time_filename = at_time_filename or ("output_at_time" in htcfile and os.path.join(htcfile.modelpath, htcfile.output_at_time.filename[0] + ".dat"))
 #         pc_filename = pc_filename or os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0])
 #         ae_filename = ae_filename or os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0])
 #         
-        mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
-        if blade_name is None:
-            blade_name = htcfile.aero.link[2]
-        self.mainbody_blade = [mb for mb in mainbodies if mb.name[0] == blade_name][0]
-        H2aeroBladeInfo.__init__(self, htcfile, ae_filename=ae_filename, pc_filename=pc_filename, at_time_filename=at_time_filename, blade_name=blade_name)
+            mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
+            if blade_name is None:
+                blade_name = htcfile.aero.link[2]
+            self.mainbody_blade = htcfile.new_htc_structure.get_subsection_by_name(blade_name)
+            st_filename = st_filename or os.path.join(htcfile.modelpath, self.mainbody_blade.timoschenko_input.filename[0])
+            MainBody.__init__(self, htcfile, blade_name)
+        elif st_filename and os.path.isfile(st_filename):
+            StFile.__init__(self, st_filename)
+        H2aeroBlade.__init__(self, htcfile, ae_filename=ae_filename, pc_filename=pc_filename, at_time_filename=at_time_filename, blade_name=blade_name)
         
 #     def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
 #         
@@ -268,14 +371,12 @@ class H2BladeInfo(H2aeroBladeInfo):
 #             htcfile = HTCFile(htcfile)
 #         
 #         blade_name = blade_name or htcfile.aero.link[2]
-        st_filename = st_filename or os.path.join(htcfile.modelpath, self.mainbody_blade.timoschenko_input.filename[0])
 #         
 #         if os.path.isfile(pc_filename) and os.path.isfile(ae_filename):
 #             PCFile.__init__(self, pc_filename, ae_filename)
 #             blade_radius = self.ae_sets[1][-1,0]
-#         
-        if os.path.isfile(st_filename):
-            StFile.__init__(self, st_filename)
+#    
+        
 #         if os.path.isfile(at_time_filename):
 #             AtTimeFile.__init__(self, at_time_filename, blade_radius)
 #             self.curved_length = self.radius_curved_ac()[-1]
@@ -287,6 +388,7 @@ class H2BladeInfo(H2aeroBladeInfo):
 #         self.c2def = np.array([v.values[1:5] for v in mainbody_blade.c2_def if v.name_ == "sec"])
 #         
 #         self.hawc2_splines_data = self.hawc2_splines()       
+    
     @property
     def c2def(self):
         if not hasattr(self, "_c2def"):
@@ -294,7 +396,7 @@ class H2BladeInfo(H2aeroBladeInfo):
         return self._c2def
    
 #         
-# class H2aeroBladeInfo(H2BladeInfo):
+# class H2aeroBlade(H2Blade):
 # 
 #     def __init__(self, at_time_filename, ae_filename, pc_filename, htc_filename):
 #         """

wetb/hawc2/mainbody.py

@@ -12,14 +12,23 @@ from wetb.hawc2.st_file import StFile
 
 
 class MainBody():
-    def __init__(self, htc_filename, modelpath, body_name):
-        self.htcfile = htcfile = HTCFile(htc_filename, modelpath)
+    def __init__(self, htcfile, body_name):
+        if isinstance(htcfile, str):
+            htcfile = HTCFile(htcfile)
+        self.htcfile = htcfile
         s = htcfile.new_htc_structure
         main_bodies = {s[k].name[0]:s[k] for k in s.keys() if s[k].name_ == "main_body"}
         self.main_body = main_bodies[body_name]
+        self.main_body = s.get_subsection_by_name(body_name)
         self.stFile = StFile(os.path.join(htcfile.modelpath, self.main_body.timoschenko_input.filename[0]))
-        self.c2def = np.array([v.values[1:5] for v in self.main_body.c2_def if v.name_ == "sec"])
+        #self.c2def = np.array([v.values[1:5] for v in self.main_body.c2_def if v.name_ == "sec"])
         self.concentrated_mass = [cm.values for cm in self.main_body if cm.name_.startswith('concentrated_mass')]
+        
+    @property
+    def c2def(self):
+        if not hasattr(self, "_c2def"):
+            self._c2def = np.array([v.values[1:5] for v in self.main_body.c2_def if v.name_ == "sec"])
+        return self._c2def
 
     def plot_xz_geometry(self, plt=None):
         if plt is None:
@@ -43,7 +52,7 @@ class MainBody():
         plt.xlabel("z")
         plt.ylabel("y")
         z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
-        plt.plot(self.c2def[:, 2], self.c2def[:, 1], label='Center line')
+        plt.plot(self.c2def[:, 2], self.c2def[:, 1], '.-', label='Center line')
         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.stFile.y_e(z), label='Elastic center')
         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.stFile.y_cg(z), label='Mass center')
         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.stFile.y_sh(z), label='Shear center')
@@ -52,106 +61,106 @@ class MainBody():
         plt.legend()
 
         
-
-
-class BladeData(object):
-    def plot_xz_geometry(self, plt):
-
-        z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
-        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]), label='Center line')
-        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) + self.pcFile.chord(z) / 2, label='Leading edge')
-        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) - self.pcFile.chord(z) / 2, label="Trailing edge")
-        x, y, z = self.hawc2_splines()
-        #plt.plot(z, x, label='Hawc2spline')
-
-    def plot_yz_geometry(self, plt):
-
-        z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
-        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]), label='Center line')
-        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.pcFile.thickness(z) / 100 * self.pcFile.chord(z) / 2, label='Suction side')
-        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) - self.pcFile.thickness(z) / 100 * self.pcFile.chord(z) / 2, label="Pressure side")
-        x, y, z = self.hawc2_splines()
-        #plt.plot(z, y, label='Hawc2spline')
-
-    def hawc2_splines(self):
-        curve_z = np.r_[0, np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))]
-        curve_z_nd = curve_z / curve_z[-1]
-
-        def akima(x, y):
-            n = len(x)
-            var = np.zeros((n + 3))
-            z = np.zeros((n))
-            co = np.zeros((n, 4))
-            for i in range(n - 1):
-                var[i + 2] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
-            var[n + 1] = 2 * var[n] - var[n - 1]
-            var[n + 2] = 2 * var[n + 1] - var[n]
-            var[1] = 2 * var[2] - var[3]
-            var[0] = 2 * var[1] - var[2]
-
-            for i in range(n):
-                wi1 = abs(var[i + 3] - var[i + 2])
-                wi = abs(var[i + 1] - var[i])
-                if (wi1 + wi) == 0:
-                    z[i] = (var[i + 2] + var[i + 1]) / 2
-                else:
-                    z[i] = (wi1 * var[i + 1] + wi * var[i + 2]) / (wi1 + wi)
-
-            for i in range(n - 1):
-                dx = x[i + 1] - x[i]
-                a = (z[i + 1] - z[i]) * dx
-                b = y[i + 1] - y[i] - z[i] * dx
-                co[i, 0] = y[i]
-                co[i, 1] = z[i]
-                co[i, 2] = (3 * var[i + 2] - 2 * z[i] - z[i + 1]) / dx
-                co[i, 3] = (z[i] + z[i + 1] - 2 * var[i + 2]) / dx ** 2
-            co[n - 1, 0] = y[n - 1]
-            co[n - 1, 1] = z[n - 1]
-            co[n - 1, 2] = 0
-            co[n - 1, 3] = 0
-            return co
-
-        def coef2spline(s, co):
-            x, y = [], []
-            for i, c in enumerate(co.tolist()[:-1]):
-                p = np.poly1d(c[::-1])
-                z = np.linspace(0, s[i + 1] - s[i ], 10)
-                x.extend(s[i] + z)
-                y.extend(p(z))
-            return y
-
-        x, y, z = [coef2spline(curve_z_nd, akima(curve_z_nd, self.c2def[:, i])) for i in range(3)]
-        return x, y, z
-
-
-
-class Blade(MainBody, BladeData):
-    def __init__(self, htc_filename, modelpath=None, blade_number=1):
-        
-        self.htcfile = htcfile = HTCFile(htc_filename, modelpath)
-        
-        blade_name = [link[2] for link in htcfile.aero if link.name_.startswith('link') and link[0]==blade_number][0]
-        MainBody.__init__(self, htc_filename, modelpath, blade_name)
-        self.pcFile = PCFile(os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0]),
-                        os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0]))
-        
-    def plot_xz_geometry(self, plt=None):
-        if plt is None:
-            import matplotlib.pyplot as plt
-            plt.figure()
-
-        MainBody.plot_xz_geometry(self, plt)
-        BladeData.plot_xz_geometry(self, plt=plt)
-        plt.legend()
-    
-    def plot_geometry_yz(self, plt=None):
-        if plt is None:
-            import matplotlib.pyplot as plt
-            plt.figure()
-
-        BladeData.plot_yz_geometry(self, plt=plt)
-        MainBody.plot_yz_geometry(self, plt)
-        
+# 
+# 
+# class BladeData(object):
+#     def plot_xz_geometry(self, plt):
+# 
+#         z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+#         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]), label='Center line')
+#         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) + self.pcFile.chord(z) / 2, label='Leading edge')
+#         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) - self.pcFile.chord(z) / 2, label="Trailing edge")
+#         x, y, z = self.hawc2_splines()
+#         #plt.plot(z, x, label='Hawc2spline')
+# 
+#     def plot_yz_geometry(self, plt):
+# 
+#         z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+#         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]), label='Center line')
+#         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.pcFile.thickness(z) / 100 * self.pcFile.chord(z) / 2, label='Suction side')
+#         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) - self.pcFile.thickness(z) / 100 * self.pcFile.chord(z) / 2, label="Pressure side")
+#         x, y, z = self.hawc2_splines()
+#         #plt.plot(z, y, label='Hawc2spline')
+# 
+#     def hawc2_splines(self):
+#         curve_z = np.r_[0, np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))]
+#         curve_z_nd = curve_z / curve_z[-1]
+# 
+#         def akima(x, y):
+#             n = len(x)
+#             var = np.zeros((n + 3))
+#             z = np.zeros((n))
+#             co = np.zeros((n, 4))
+#             for i in range(n - 1):
+#                 var[i + 2] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
+#             var[n + 1] = 2 * var[n] - var[n - 1]
+#             var[n + 2] = 2 * var[n + 1] - var[n]
+#             var[1] = 2 * var[2] - var[3]
+#             var[0] = 2 * var[1] - var[2]
+# 
+#             for i in range(n):
+#                 wi1 = abs(var[i + 3] - var[i + 2])
+#                 wi = abs(var[i + 1] - var[i])
+#                 if (wi1 + wi) == 0:
+#                     z[i] = (var[i + 2] + var[i + 1]) / 2
+#                 else:
+#                     z[i] = (wi1 * var[i + 1] + wi * var[i + 2]) / (wi1 + wi)
+# 
+#             for i in range(n - 1):
+#                 dx = x[i + 1] - x[i]
+#                 a = (z[i + 1] - z[i]) * dx
+#                 b = y[i + 1] - y[i] - z[i] * dx
+#                 co[i, 0] = y[i]
+#                 co[i, 1] = z[i]
+#                 co[i, 2] = (3 * var[i + 2] - 2 * z[i] - z[i + 1]) / dx
+#                 co[i, 3] = (z[i] + z[i + 1] - 2 * var[i + 2]) / dx ** 2
+#             co[n - 1, 0] = y[n - 1]
+#             co[n - 1, 1] = z[n - 1]
+#             co[n - 1, 2] = 0
+#             co[n - 1, 3] = 0
+#             return co
+# 
+#         def coef2spline(s, co):
+#             x, y = [], []
+#             for i, c in enumerate(co.tolist()[:-1]):
+#                 p = np.poly1d(c[::-1])
+#                 z = np.linspace(0, s[i + 1] - s[i ], 10)
+#                 x.extend(s[i] + z)
+#                 y.extend(p(z))
+#             return y
+# 
+#         x, y, z = [coef2spline(curve_z_nd, akima(curve_z_nd, self.c2def[:, i])) for i in range(3)]
+#         return x, y, z
+
+
+
+# class Blade(MainBody, BladeData):
+#     def __init__(self, htc_filename, modelpath=None, blade_number=1):
+#         
+#         self.htcfile = htcfile = HTCFile(htc_filename, modelpath)
+#         
+#         blade_name = [link[2] for link in htcfile.aero if link.name_.startswith('link') and link[0]==blade_number][0]
+#         MainBody.__init__(self, htc_filename, modelpath, blade_name)
+#         self.pcFile = PCFile(os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0]),
+#                         os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0]))
+#         
+#     def plot_xz_geometry(self, plt=None):
+#         if plt is None:
+#             import matplotlib.pyplot as plt
+#             plt.figure()
+# 
+#         MainBody.plot_xz_geometry(self, plt)
+#         BladeData.plot_xz_geometry(self, plt=plt)
+#         plt.legend()
+#     
+#     def plot_geometry_yz(self, plt=None):
+#         if plt is None:
+#             import matplotlib.pyplot as plt
+#             plt.figure()
+# 
+#         BladeData.plot_yz_geometry(self, plt=plt)
+#         MainBody.plot_yz_geometry(self, plt)
+#         
 
 
 

wetb/hawc2/pc_file.py

@@ -16,7 +16,7 @@ standard_library.install_aliases()
 from wetb.hawc2.ae_file import AEFile
 import numpy as np
 
-class PCFile(AEFile):
+class PCFile(object):
     """Read HAWC2 PC (profile coefficients) file
 
     examples
@@ -37,9 +37,7 @@ class PCFile(AEFile):
     >>> pcfile.CM(36,10) # CM at radius=36m and AOA=10deg
     -0.1103
     """
-    def __init__(self, filename, ae_filename):
-        AEFile.__init__(self, ae_filename)
-
+    def __init__(self, filename):
         with open (filename) as fid:
             lines = fid.readlines()
         nsets = int(lines[0].split()[0])
@@ -61,9 +59,7 @@ class PCFile(AEFile):
                 lptr += n_rows
             self.pc_sets[nset] = (np.array(thicknesses), profiles)
 
-    def _Cxxx(self, radius, alpha, column, ae_set_nr=1):
-        thickness = self.thickness(radius, ae_set_nr)
-        pc_set_nr = self.pc_set_nr(radius, ae_set_nr)
+    def _Cxxx(self, thickness, alpha, column, ae_set_nr=1, pc_set_nr=1):
         thicknesses, profiles = self.pc_sets[pc_set_nr]
         index = np.searchsorted(thicknesses, thickness)
         if index == 0:
@@ -75,9 +71,7 @@ class PCFile(AEFile):
         th0, th1 = thicknesses[index - 1:index + 1]
         return Cx0 + (Cx1 - Cx0) * (thickness - th0) / (th1 - th0)
     
-    def _CxxxH2(self, radius, alpha, column, ae_set_nr=1):
-        thickness = self.thickness(radius, ae_set_nr)
-        pc_set_nr = self.pc_set_nr(radius, ae_set_nr)
+    def _CxxxH2(self, thickness, alpha, column, ae_set_nr=1,pc_set_nr=1):
         thicknesses, profiles = self.pc_sets[pc_set_nr]
         index = np.searchsorted(thicknesses, thickness)
         if index == 0:
@@ -95,13 +89,13 @@ class PCFile(AEFile):
     
         
 
-    def CL(self, radius, alpha, ae_set_nr=1):
+    def CL(self, thickness, alpha, ae_set_nr=1,pc_set_nr=1):
         """Lift coefficient
 
         Parameters
         ---------
-        radius : float
-            radius [m]
+        thickness : float
+            thickness [5]
         alpha : float
             Angle of attack [deg]
         ae_set_nr : int optional
@@ -111,13 +105,13 @@ class PCFile(AEFile):
         -------
         Lift coefficient : float
         """
-        return self._Cxxx(radius, alpha, 1, ae_set_nr)
+        return self._Cxxx(thickness, alpha, 1, ae_set_nr, pc_set_nr)
 
 
-    def CL_H2(self, radius, alpha, ae_set_nr=1):
-        return self._CxxxH2(radius, alpha, 1, ae_set_nr)
+    def CL_H2(self, thickness, alpha, ae_set_nr=1, pc_set_nr=1):
+        return self._CxxxH2(thickness, alpha, 1, ae_set_nr, pc_set_nr)
     
-    def CD(self, radius, alpha, ae_set_nr=1):
+    def CD(self, thickness, alpha, ae_set_nr=1, pc_set_nr=1):
         """Drag coefficient
 
         Parameters
@@ -133,10 +127,10 @@ class PCFile(AEFile):
         -------
         Drag coefficient : float
         """
-        return self._Cxxx(radius, alpha, 2, ae_set_nr)
+        return self._Cxxx(thickness, alpha, 2, ae_set_nr, pc_set_nr)
 
-    def CM(self, radius, alpha, ae_set_nr=1):
-        return self._Cxxx(radius, alpha, 3, ae_set_nr)
+    def CM(self, thickness, alpha, ae_set_nr=1, pc_set_nr=1):
+        return self._Cxxx(thickness, alpha, 3, ae_set_nr,pc_set_nr)
 
 if __name__ == "__main__":
     pc = PCFile(r"C:\mmpe\Projects\inflow\Hawc2aero_setup/data/Hawc_pc.b52", r"C:\mmpe\Projects\inflow\Hawc2aero_setup/data/S36_ae_h2.001")

wetb/hawc2/tests/test_blade.py

+'''
+Created on 3. maj 2017
+
+@author: mmpe
+'''
+import os
+import unittest
+
+from wetb.hawc2.blade import H2Blade
+import numpy as np
+
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/')  # test file path
+class Test(unittest.TestCase):
+
+
+#     def testBladeInfo(self):
+#         bi = H2Blade(tfp + "simulation_setup/DTU10MWRef6.0/htc/DTU_10MW_RWT.htc")
+#         if 0:
+#             import matplotlib.pyplot as plt
+#             print (dir(bi))
+#             #print (bi.radius_s())
+#             plt.plot(bi.radius_s(), bi.twist())
+#             plt.plot(bi.c2def[:,2], bi.c2def[:,3])
+#             x = np.linspace(0,1,1000)
+#             plt.plot(bi.blade_radius*x, bi.c2nd(x)[:,3])
+#             plt.show()
+
+            
+    def testBladeInfo_AE(self):
+        bi = H2Blade(None, tfp + "NREL_5MW_ae.txt", tfp + "NREL_5MW_pc.txt")
+        self.assertEqual(bi.thickness(32), 23.78048780487805)
+        self.assertEqual(bi.chord(32), 3.673)
+        self.assertEqual(bi.pc_set_nr(32), 1)
+
+    def testBladeInfo_PC_AE(self):
+        bi = H2Blade(None, tfp + "NREL_5MW_ae.txt", tfp + "NREL_5MW_pc.txt")
+        self.assertEqual(bi.CL(36, 10), 1.358)
+        self.assertEqual(bi.CD(36, 10), 0.0255)
+        self.assertEqual(bi.CM(36, 10), -0.1103)
+            
+#     def test_curved_length2radius(self):
+#         bi = H2BladeInfo(tfp + "simulation_setup/DTU10MWRef6.0/htc/DTU_10MW_RWT.htc")
+#         np.testing.assert_array_equal(bi.xyztwist(86.4979, curved_length=True), bi.xyztwist(86.3655))
+        
+            
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testBladeInfo']
+    unittest.main()
\ No newline at end of file

wetb/hawc2/tests/test_main_body.py

@@ -5,13 +5,14 @@ Created on 01/08/2016
 '''
 import unittest
 import os
-from wetb.hawc2.mainbody import MainBody, Blade
+from wetb.hawc2.mainbody import MainBody
+from wetb.hawc2.blade import H2Blade
 
 tfp = os.path.join(os.path.dirname(__file__), 'test_files/')  # test file path
 class TestMainBody(unittest.TestCase):
     
     def test_MainBody(self):
-        mainbody = MainBody(tfp+"htcfiles/test.htc", "../", "towertop")
+        mainbody = MainBody(tfp+"htcfiles/test.htc", "blade1")
         if 0:
             import matplotlib.pyplot as plt
             plt.figure()
@@ -21,7 +22,7 @@ class TestMainBody(unittest.TestCase):
             plt.show()
              
     def test_Blade(self):
-        blade = Blade(tfp+"htcfiles/test.htc", "../")
+        blade = H2Blade(tfp+"htcfiles/test.htc")
         if 0:
             import matplotlib.pyplot as plt
             plt.figure()

wetb/hawc2/tests/test_pc_file.py

@@ -8,6 +8,7 @@ from __future__ import print_function
 from __future__ import division
 from __future__ import absolute_import
 from future import standard_library
+from wetb.hawc2.ae_file import AEFile
 standard_library.install_aliases()
 import os
 import unittest
@@ -24,19 +25,13 @@ class TestPCFile(unittest.TestCase):
         unittest.TestCase.setUp(self)
         self.testfilepath = os.path.join(os.path.dirname(__file__), 'test_files/')  # test file path
 
-
     def test_PCFile_ae(self):
-        pc = PCFile(self.testfilepath + "NREL_5MW_pc.txt", self.testfilepath + "NREL_5MW_ae.txt")
-        self.assertEqual(pc.thickness(32), 23.78048780487805)
-        self.assertEqual(pc.chord(32), 3.673)
-        self.assertEqual(pc.pc_set_nr(32), 1)
-
-    def test_PCFile2(self):
-        pc = PCFile(self.testfilepath + "NREL_5MW_pc.txt", self.testfilepath + "NREL_5MW_ae.txt")
-        self.assertEqual(pc.CL(36, 10), 1.358)
-        self.assertEqual(pc.CD(36, 10), 0.0255)
-        self.assertEqual(pc.CM(36, 10), -0.1103)
-
+        pc = PCFile(self.testfilepath + "NREL_5MW_pc.txt")
+        ae = AEFile(self.testfilepath + "NREL_5MW_ae.txt")
+        thickness = ae.thickness(36)
+        self.assertEqual(pc.CL(thickness, 10), 1.358)
+        self.assertEqual(pc.CD(thickness, 10), 0.0255)
+        self.assertEqual(pc.CM(thickness, 10), -0.1103)
 
 if __name__ == "__main__":
     #import sys;sys.argv = ['', 'Test.testName']