new pc file without ae file.

Old pc-ae file integrated in bladeinfo

wetb/hawc2/bladeinfo.py

@@ -11,6 +11,7 @@ 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
 
 
 
@@ -61,7 +62,7 @@ from wetb.hawc2.st_file import StFile
 #         """
 
 
-class H2aeroBladeInfo(PCFile, AtTimeFile):
+class H2aeroBladeInfo(PCFile, AEFile, AtTimeFile):
     """Provide HAWC2 info about a blade
     
     From AE file:
@@ -86,36 +87,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,6 +224,83 @@ 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):
     """Provide HAWC2 info about a blade
@@ -247,18 +326,20 @@ 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]
+            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])
         H2aeroBladeInfo.__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,13 +349,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):
+        if st_filename and 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)

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_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']