ktreff_aerofoil.py

import math

import matlab.engine
import numpy as np
import pandas as pd
from scipy.interpolate import griddata


class KarmanTrefftzAerofoil(object):
    """
    Object to represent Karman-Trefftz aerofoil generated by MATLAB.
    """
    def __init__(self, **kwargs):
        self.alpha = float(kwargs['alpha']) * math.pi / 180.0   # Input in degrees, MATLAB requires radians
        self.eps = float(kwargs['eps'])
        self.beta = float(kwargs['beta'])
        self.tau = float(kwargs['tau'])
        self.n = int(kwargs['n'])
        self.data = np.empty((self.n, 3))
        self.Cl = 0.0
        self.lower = None
        self.upper = None
        self.camber = None
        self.thickness = None
        self.dpan = True if 'dpan' in kwargs.keys() else False

    def __call__(self, *args, **kwargs):
        self._get_ktreff_data_matlab()
        self._calc_lower_upper()
        self._calc_camber()
        self._calc_thickness()
        if self.dpan:
            self._get_doublet_data_matlab()

    def _get_ktreff_data_matlab(self):
        """
        Function to call MATLAB ktreff function and store data in object.
        """
        eng = matlab.engine.start_matlab()
        x, y, cp, cl = eng.ktreff(self.alpha, self.eps, self.beta, self.tau, float(self.n - 1), nargout=4)
        self.data[:, 0] = x[0]
        self.data[:, 1] = y[0]
        self.data[:, 2] = cp[0]
        self.Cl = cl

    def _get_doublet_data_matlab(self):
        """
        Function to call MATLAB dpan function and store data in object.
        """
        x = matlab.double(self.data[:, 0].tolist())[0]
        y = matlab.double(self.data[:, 1].tolist())[0]
        eng = matlab.engine.start_matlab()
        cp, cl = eng.dpan(float(self.n - 1), self.alpha, x, y, nargout=2)
        self.dpan_cl = cl
        self.dpan_cp = cp[0]

    def get_excel(self):
        """
        Function to write .xlsx file containing x, y, camber, thickness and Cp.
        """
        data = pd.DataFrame(columns=['x', 'y', 'Cp'])
        data['x'] = self.data[:, 0]
        data['y'] = self.data[:, 1]
        data['Cp'] = self.data[:, 2]

        camber = pd.DataFrame(columns=['x', 'camber'])
        camber['x'] = self.camber[:, 0]
        camber['camber'] = self.camber[:, 1]

        thickness = pd.DataFrame(columns=['thickness'])
        thickness['thickness'] = self.thickness
        return data, camber, thickness

    def write_dat(self):
        """
        Function to write .dat file containing x, y coordinates.
        """
        data = self.get_data(True)
        data.savetxt('{}.dat'.format(self.get_name()))

    def get_data(self, coords_only):
        """
        :param coords_only: True if only coordinate data required.
        :return: Return array of data for aerofoil.
        """
        return self.data[:, :-1] if coords_only else self.data

    def _calc_lower_upper(self):
        data = self.get_data(True)
        lower = data[0:self.n/2]
        upper = data[self.n/2:self.n]
        self.lower = lower[lower[:, 0].argsort()][::-1]
        self.upper = upper[upper[:, 0].argsort()][::-1]

    def get_lower_upper(self):
        return self.lower, self.upper

    def _calc_camber(self):
        self.camber = 0.5 * (self.lower + self.upper)

    def get_camber(self):
        return self.camber

    def get_max_camber(self):
        return self.camber[np.argmax(self.camber[:, 1])]

    def _calc_thickness(self):
        interp_upper = griddata(self.upper[:, 0], self.upper[:, 1], self.lower[:, 0], method='linear')
        thickness = interp_upper - self.lower[:, 1]
        thickness = thickness[np.where(np.logical_not(np.isnan(thickness)))]
        self.thickness = thickness

    def get_thickness(self):
        return self.thickness

    def get_max_thickness(self):
        return np.max(self.thickness)

    def get_naca_digits(self):
        max_camber_x, max_camber = self.get_max_camber()
        max_thickness = self.get_max_thickness()
        return int(round(max_camber * 100)), int(round(max_camber_x * 10)), int(round(max_thickness * 100))

    def get_name(self):
        return 'ktreff(alpha={},eps={},beta={},tau={},n={})'.format(np.round(self.alpha * math.pi / 180.0, 1),
                                                                    self.eps, self.beta, self.tau,
                                                                    self.n)