blade.py

'''
Created on 3. maj 2017

@author: mmpe
'''
import numpy as np
from wetb.hawc2.at_time_file import AtTimeFile
from wetb.hawc2.pc_file import PCFile
from wetb.hawc2.htc_file import HTCFile
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



# class BladeInfo(object):
#     def twist(self, radius=None):
#         """Aerodynamic twist [deg]. Negative when leading edge is twisted towards wind(opposite to normal definition)\n
# 
#         Parameters
#         ---------
#         radius : float or array_like, optional
#             radius [m] of interest
#             If None (default) the twist of all points are returned
# 
#         Returns
#         -------
#         twist [deg] : float
#         """
#         raise NotImplementedError()
# 
#     def chord(self, radius=None):
#         """Aerodynamic chord
# 
#         Parameters
#         ---------
#         radius : float or array_like, optional
#             radius [m] of interest
#             If None (default) the twist of all points are returned
# 
#         Returns
#         -------
#         chord [m] : float
#         """
#         raise NotImplementedError()
# 
#     
# 
#     def c2nd(self, radius_nd):
#         """Return center line position
# 
#         Parameters
#         ---------
#         radius_nd : float or array_like, optional
#             non dimensional radius
# 
#         Returns
#         -------
#         x,y,z : float
#         """


class H2aeroBlade(PCFile, AEFile, AtTimeFile):
    """Provide HAWC2 info about a blade
    
    From AE file:
    - chord(radius=None, set_nr=1):
    - thickness(radius=None, set_nr=1)
    - radius_ae(radius=None, set_nr=1)
    
    From PC file
    - CL(radius, alpha, ae_set_nr=1)
    - CD(radius, alpha, ae_set_nr=1)
    - CM(radius, alpha, ae_set_nr=1)
    
    From at_time_filename
    - attribute_names
    - xxx(radius=None, curved_length=None) # xxx for each attribute name
    - radius_curved_ac(radius=None) # Curved length of nearest/all aerodynamic calculation points
    
    From ST file
    - radius_st(radius=None, mset=1, set=1)
    - xxx(radius=None, mset=1, set=1) # xxx for each of r, m, x_cg,y_cg, ri_x, ri_y, xs, ys, E, G, Ix, Iy, K, kx, ky, A, pitch, xe, ye
    with template
        
    """

    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):
            AEFile.__init__(self, ae_filename)
            PCFile.__init__(self, pc_filename)
            blade_radius = self.ae_sets[1][-1,0]
        
        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:
            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]

        
        
        
        
        
    @property
    def c2def(self):
        if not hasattr(self, "_c2def"):
            self._c2def = np.array([v.values[1:5] for v in self.htcfile.blade_c2_def if v.name_ == "sec"])
        return self._c2def
    
        
    def c2nd(self, l_nd, curved_length=True):
        curve_l_nd, x, y, z, twist = self.hawc2_splines_data
        if curved_length:
            l_nd = (np.max([np.min([l_nd, np.ones_like(l_nd)], 0), np.zeros_like(l_nd) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
            return np.array([np.interp(l_nd, curve_l_nd, xyz) for xyz in [x,y,z, twist]]).T
        else:
            l_nd = (np.max([np.min([l_nd, np.ones_like(l_nd)], 0), np.zeros_like(l_nd)], 0))
            return np.array([np.interp(l_nd, z/z[-1], xyz) for xyz in [x,y,z, twist]]).T
            
    def c2(self, l, curved_length=True):
        if curved_length:
            L = self.curved_length
        else:
            L = self.c2def[-1,2]
        return self.c2nd(l/L, curved_length)
    
    def hawc2_splines(self):
        curve_l = np.r_[0, np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))]
        curve_l_nd = curve_l / curve_l[-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, endpoint=i >= co.shape[0] - 2)
                x.extend(s[i] + z)
                y.extend(p(z))
            return x,y

        x, y, z, twist = [coef2spline(curve_l_nd, akima(curve_l_nd, self.c2def[:, i]))[1] for i in range(4)]
        curve_l_nd = coef2spline(curve_l_nd, akima(curve_l_nd, self.c2def[:, 0]))[0]
        return curve_l_nd, x, y, z, twist

    def xyztwist(self, l=None, curved_length=False):
        """Return splined x,y,z and twist 
        
        Parameters
        ----------
        l : int, float, arraylike or None, optional
            Position of interest, seee curved_length\n
            If None (default) all x, y, z, and twist defined in c2def
        curved_length : bool, optional
            - If False: l is z coordinate of section
            - If True: l is curved length
            
        Returns
        -------
        x,y,z,twist
                """
        if l is None:
            return self.c2def[:, :3]
        else:
            r_nd = np.linspace(0,1,1000)
            if curved_length:
                #curved_length = np.cumsum(np.sqrt((np.diff(self.c2nd(np.linspace(0,1,100)),1,0)[:,:3]**2).sum(1)))
                return self.c2nd(l/self.radius_curved_ac()[-1])    
#                 z_nd = (np.cos(np.linspace(np.pi, np.pi*2,len(curved_length)-1))+1)/2
#                 assert np.all(l>=curved_length[0]) and np.all(l<=curved_length[-1])
#                 return self.c2nd(r_nd[np.argmin(np.abs(curved_length-l))+1])    
            else:
                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)
            
            
    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:
    - chord(radius=None, set_nr=1):
    - thickness(radius=None, set_nr=1)
    - radius_ae(radius=None, set_nr=1)
    
    From PC file
    - CL(radius, alpha, ae_set_nr=1)
    - CD(radius, alpha, ae_set_nr=1)
    - CM(radius, alpha, ae_set_nr=1)
    
    From at_time_filename
    - attribute_names
    - xxx(radius=None, curved_length=None) # xxx for each attribute name
    - radius_curved_ac(radius=None) # Curved length of nearest/all aerodynamic calculation points
    
    From ST file
    - radius_st(radius=None, mset=1, set=1)
    - xxx(radius=None, mset=1, set=1) # xxx for each of r, m, x_cg,y_cg, ri_x, ri_y, xs, ys, E, G, Ix, Iy, K, kx, ky, A, pitch, xe, ye
    with template
        
    """
    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 = 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):
#         
#         if isinstance(htcfile, str):
#             assert htcfile.lower().endswith('.htc')
#             htcfile = HTCFile(htcfile)
#         
#         blade_name = blade_name or htcfile.aero.link[2]
#         
#         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(at_time_filename):
#             AtTimeFile.__init__(self, at_time_filename, blade_radius)
#             self.curved_length = self.radius_curved_ac()[-1]
#         else:
#             raise NotImplementedError
#             #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.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"):
            self._c2def = np.array([v.values[1:5] for v in self.mainbody_blade.c2_def if v.name_ == "sec"])
        return self._c2def
   
#         
# class H2aeroBlade(H2Blade):
# 
#     def __init__(self, at_time_filename, ae_filename, pc_filename, htc_filename):
#         """
#         at_time_filename: file name of at time file containing twist and chord data
#         """
#         PCFile.__init__(self, pc_filename, ae_filename)
#         blade_radius = self.ae_sets[1][-1,0]
#         AtTimeFile.__init__(self, at_time_filename, blade_radius)
# 
#         assert('twist' in self.attribute_names)
#         htcfile = HTCFile(htc_filename)
# 
# 
#         self.c2def = np.array([v.values[1:5] for v in htcfile.blade_c2_def if v.name_ == "sec"])
#         #self._c2nd = interp1d(self.c2def[:, 2] / self.c2def[-1, 2], self.c2def[:, :3], axis=0, kind='cubic')
# 
#         ac_radii = self.radius_curved_ac()
#         self.hawc2_splines_data = self.hawc2_splines()
#         self.curved_length = np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))[-1]
#         
# #         c2_axis_length = np.r_[0, np.cumsum(np.sqrt((np.diff(self.c2def[:, :3], 1, 0) ** 2).sum(1)))]
# #         self._c2nd = interp1d(c2_axis_length / c2_axis_length[-1], self.c2def[:, :3], axis=0, kind='cubic')
#         #self._c2nd = interp1d(self.c2def[:,2]/self.c2def[-1,2], self.c2def[:, :3], axis=0, kind='cubic')
# # 
# #     def c2nd(self, r_nd):
# #         r_nd_min = np.zeros_like(r_nd) + self.c2def[0, 2] / self.c2def[-1, 2]
# #         r_nd_max = np.ones_like(r_nd)
# #         r_nd = np.max([np.min([r_nd, r_nd_max], 0), r_nd_min], 0)
# #         return self._c2nd(r_nd)
#