diff --git a/.gitignore b/.gitignore
index 409444e57b1137c057d38885f6e3ff7c3d01b282..62030f021ef1e5e8be30d24975d1b2be3d27d0b6 100644
--- a/.gitignore
+++ b/.gitignore
@@ -86,6 +86,7 @@ nosetests.xml
 
 *.xyz
 *.xyz.fvbnd
+*.x2d
 
 Icon?
 ehthumbs.db
diff --git a/PGL/components/airfoil.py b/PGL/components/airfoil.py
index 899a31c4270623aead50a85f1d29ce6b7cf925bb..cb6d7fc0389d23bf4d02ab9cd60d5de71ee6903b 100644
--- a/PGL/components/airfoil.py
+++ b/PGL/components/airfoil.py
@@ -8,8 +8,8 @@ from scipy.optimize import minimize
 from scipy.interpolate import pchip, Akima1DInterpolator, interp1d
 
 from PGL.main.geom_tools import curvature, calculate_length
-from PGL.main.curve import Curve
-from PGL.main.bezier import FitBezier
+from PGL.main.curve import Curve, Line
+from PGL.main.bezier import FitBezier, BezierCurve
 from PGL.main.naturalcubicspline import NaturalCubicSpline
 
 
@@ -111,44 +111,76 @@ class AirfoilShape(Curve):
             optional trailing edge cell size. If set to -1 the cell size will increase
             from the LE to TE according to the tanh distribution function used
             in distfunc
-        close_te: int
-            close trailing edge with line segment consisting of n points.
+        close_te: int/bool
+            close trailing edge with line segment consisting of n points if an integer,
+            if bool and True the number of points will be determined automatically, such the
+            grid spacing is similar to the one at the LE 
+            
         """
-
+        # trailing edge length
+        self.main_seg = Curve(points=self.points, spline=self.spline_type, compute_dp=False)
+        len_TE = np.sqrt(((self.points[-1,:] - self.points[0,:])**2).sum()) / self.main_seg.smax
+        # check whether TE is already closed
+        if len_TE < 1e-10:
+            print('TE is closed already and cannot be closed!')
+            close_te = 0
+        # automatic determination of cells to be used to close TE, such that it is similar to the 
+        # one used at the LE, if close_te=False it is not closed 
+        if isinstance(close_te, (bool)):
+            if close_te:
+                if dLE:
+                    ds_LE = dLE
+                else:
+                    ds_LE = self.leading_edge_dist(ni)
+                # estimate number of points on TE, such that the grid spacing is close to LE 
+                ni_TE = int(np.round(len_TE/ds_LE))
+                # ensure there to be one point at least 
+                if ni_TE ==0:
+                    ni_TE = 1
+                close_te = ni_TE    
+            else:
+                close_te = 0
+        # close trailing edge
         if close_te > 0:
             if close_te % 2 == 0:
                 close_te += 1
 
-            # construct the trailing edge segments
-            lte_seg = np.array([np.linspace(self.TE.copy()[i], self.points[0, i], close_te//2+1) for i in range(self.nd)]).T
-            ute_seg = np.array([np.linspace(self.points[-1, i], self.TE.copy()[i], close_te//2+1) for i in range(self.nd)]).T
-            lte_seg = Curve(points=lte_seg, spline='linear')
-            ute_seg = Curve(points=ute_seg, spline='linear')
+            # construct the trailing edge segment, from suction to pressure side
+            self.ni_TE = close_te
+            self.te_seg = Line(self.main_seg.points[-1,:], self.main_seg.points[0,:], max(2, self.ni_TE))
+            self.te_isplit = self.ni_TE//2      
 
-            main_seg = Curve(points=self.points, spline=self.spline_type, compute_dp=False)
-            minTE = lte_seg.ds[-1]/ main_seg.smax
+            minTE = self.te_seg.ds[-1]/ self.main_seg.smax
             main_ni = ni - close_te + 1
             if dLE:
-                dist_LE = dLE
+                ds_LE = dLE
             else:
-                dist_LE = self.leading_edge_dist(main_ni)
+                ds_LE = self.leading_edge_dist(main_ni)
             # include trailing/leading edge control points
             disttot = [[0., minTE, 1],
-                    [self.sLE, dist_LE, (main_ni-1)// 2 + 1],
+                    [self.sLE, ds_LE, (main_ni-1)// 2 + 1],
                     [1., minTE, main_ni]]
-            # include additional uder-defined control points and re-order
+            # include additional under-defined control points and re-order
             if(dist!=None):
                 for elem in dist:
                     disttot.append(elem)
                 disttot.sort(key=lambda x: float(x[0]))
             # apply total distribution
-            main_seg.redistribute(dist=disttot, linear=linear)
-            points = lte_seg.points.copy()
-            points = np.append(points, main_seg.points[1:], axis=0)
-            points = np.append(points, ute_seg.points[1:], axis=0)
-            # self.initialize(points)
+            self.main_seg.redistribute(dist=disttot, linear=linear)
+            # Create the TE closure
+            points = np.append(self.te_seg.points[self.te_isplit:], self.main_seg.points[1:], axis=0)
+            points = np.append(points, self.te_seg.points[1:self.te_isplit+1], axis=0)
+            # Outputs
             self.points = points
             self.ni = ni
+            # some additional outputs needed for LER generation
+            self.main_ni = main_ni
+            self.ds_TE = len_TE / self.ni_TE * self.main_seg.smax
+            self.ds_LE = ds_LE * self.main_seg.smax
+            # ensure the curve location of the LE is matching between the main segment and the profile
+            # (a very small difference would otherwise persist) 
+            self.main_seg.sLE = self.sLE
+            #self.initialize(points)
         else:
             if dist == None:
                 if even:
@@ -363,6 +395,973 @@ class AirfoilShape(Curve):
             ix = np.where(np.diff(s)>0)[0]
             spl = NaturalCubicSpline(s[ix], y[ix])
             return spl(x)
+            
+    # Implementation of different LER types
+    def bite(self, bite):
+        """
+        Create over- (bite > 0) and underbite (bite < 0)
+
+        bite: float [m]
+            size of over or underbite
+        """
+        if (bite > 0.):
+            self.overbite(bite)
+        else:
+            self.underbite(abs(bite))
+
+    def overbite(self, overbite):
+        """
+        Create overbite by moving the pressure side downstream
+
+        overbite: float [m]
+            size of overbite
+        """
+        iLE = np.argmin(self.main_seg.points[:,0])
+        points_lower = self.main_seg.points[:iLE+1,:].copy()
+        points_upper = self.main_seg.points[iLE:,:].copy()
+        # Determine the no. of points needed on the overbite
+        ni_over = max(4,int(np.ceil(abs(overbite)/self.ds_LE)/2))+2
+        if ni_over%2 != 0:
+            ni_over += 1
+        ni_upper = points_lower.shape[0]-ni_over/2
+        ni_lower = points_upper.shape[0]-ni_over/2
+        ni_over += 1
+        self.ni_over = ni_over
+        # Move lower part of airfoil streamwise
+        points_lower[:,0] += overbite
+        low = Curve(points=points_lower, spline=self.spline_type, compute_dp=False)
+        up = Curve(points=points_upper, spline=self.spline_type, compute_dp=False)
+        # Find the TE of the shifted airfoil
+        yorig_TE = self.main_seg.points[0,1]
+        xorig_TE = self.main_seg.points[0,0]
+        yrot_TE = 1.
+        ylim = 1e-5
+        dth = ylim/100.*180./np.pi
+        th = 0.
+
+        # Rotate the spline about the LE
+        while abs(yorig_TE - yrot_TE)>ylim:
+            low.rotate_z(-th,center=points_lower[-1,:])
+            self.fit = self._splcls(low.points[::-1, 0], low.points[::-1, 1])
+            yrot_TE = self.fit(xorig_TE)
+            th += dth
+        # Find the s location for that x position
+        fit = self._splcls(low.points[::-1, 0], low.s[::-1])
+        srot_TE = fit(xorig_TE)
+
+        # Change the point distribution to fit the new curve length
+        # The first part is only a dummy section which will be removed later
+        ni_dum = 10
+
+        ds = min(self.ds_LE,abs(overbite)/ni_over)
+
+        dist = [[0., self.ds_TE/low.smax, 1],
+                [srot_TE, self.ds_TE/low.smax, 1+ni_dum],
+                [1., ds/low.smax, ni_lower+ni_dum]]
+        low.redistribute(dist=dist)
+        dist = [[0, ds/up.smax, 1],
+                [1., self.ds_TE/up.smax, ni_upper]]
+        up.redistribute(dist=dist)
+        # Create trailing edge
+        self.te_seg = Line(up.points[-1,:], low.points[ni_dum,:], self.ni_TE+2)
+
+        # Create overhang
+        self.over_seg = Line(low.points[-1,:], up.points[0,:], ni_dum)
+        self.over_seg.redistribute(dist=[[0., ds/self.over_seg.smax, 1],[1., ds/self.over_seg.smax, ni_over]])
+
+        self.points_lower = points_lower
+        self.points_lower_rot = low.points
+        self.low = low
+        self.points_upper = up.points
+
+        # Assemble mdofied profile
+        points = np.append(self.te_seg.points[self.te_isplit:],low.points[1+ni_dum:], axis=0)
+        points = np.append(points,self.over_seg.points[1:], axis=0)
+        points = np.append(points,up.points[1:], axis=0)
+        points = np.append(points,self.te_seg.points[1:self.te_isplit+1], axis=0)
+        self.points = points
+        self.initialize(self.points)
+
+        return self
+
+    def underbite(self, underbite):
+        """
+        Create underbite by moving the suction side downstream
+
+        underbite: float [m]
+            size of underbite
+        """
+        iLE = np.argmin(self.main_seg.points[:,0])
+        points_lower = self.main_seg.points[:iLE+1,:].copy()
+        points_upper = self.main_seg.points[iLE:,:].copy()
+        # Determine the no. of points needed on the overbite
+        ni_over = max(4,int(np.ceil(abs(underbite)/self.ds_LE)/2))+2
+        if ni_over%2 != 0:
+            ni_over += 1
+        ni_upper = int(points_lower.shape[0]-ni_over/2)
+        ni_lower = int(points_upper.shape[0]-ni_over/2)
+        ni_over += 1
+        self.ni_over = ni_over
+        # Move lower part of airfoil streamwise
+        points_upper[:,0] += underbite
+        low = Curve(points=points_lower, spline=self.spline_type, compute_dp=False)
+        up = Curve(points=points_upper, spline=self.spline_type, compute_dp=False)
+        self.low = low
+        self.up = up
+
+        # Find the TE of the shifted airfoil
+        yorig_TE = self.main_seg.points[-1,1]
+        xorig_TE = self.main_seg.points[-1,0]
+        yrot_TE = 1.
+        ylim = 1e-5
+        dth = ylim/100.*180./np.pi
+        th = 0.
+        # Rotate the spline about the LE
+        while abs(yorig_TE - yrot_TE)>ylim:
+            up.rotate_z(-th,center=points_upper[0,:])
+            self.fit = self._splcls(up.points[:, 0], up.points[:, 1])
+            yrot_TE = self.fit(xorig_TE)
+            th += dth
+        # Find the s location for that x position
+        fit = self._splcls(up.points[:, 0], up.s)
+        srot_TE = fit(xorig_TE)
+        self.up2 = up
+
+        # Change the point distribution to fit the new curve length
+        # The first part is only a dummy section which will be removed later
+        ni_dum = 10
+
+        ds = min(self.ds_LE,abs(underbite)/ni_over)
+
+        dist = [[0, self.ds_TE/low.smax, 1],
+                [1., ds/low.smax, ni_lower]]
+        low.redistribute(dist=dist)
+        dist = [[0, ds/up.smax, 1],
+                [srot_TE, self.ds_TE/up.smax, ni_upper],
+                [1., self.ds_TE/up.smax, ni_upper+ni_dum]]
+        up.redistribute(dist=dist)
+
+        # Create trailing edge
+        self.te_seg = Line(up.points[ni_upper-1,:], low.points[0,:], self.ni_TE+2)
+
+        # Create overhang
+        self.over_seg = Line(low.points[-1,:], up.points[0,:], ni_dum)
+
+        self.over_seg.redistribute(dist=[[0., ds/self.over_seg.smax, 1],[1., ds/self.over_seg.smax, ni_over]])
+
+        self.points_upper = points_upper
+        self.points_upper_rot = up.points
+        self.low = low
+        self.points_lower = low.points
+        self.ni_upper = ni_upper
+
+        # Assemble mdofied profile
+        points = np.append(self.te_seg.points[self.te_isplit:],low.points[1:], axis=0)
+        points = np.append(points,self.over_seg.points[1:], axis=0)
+        points = np.append(points,up.points[1:ni_upper], axis=0)
+        points = np.append(points,self.te_seg.points[1:self.te_isplit+1], axis=0)
+        self.points = points
+        self.initialize(self.points)
+
+        return self
+
+    def wavy(self, s_start, s_end, act_len, amp, lam, ni):
+        """
+        Create sinusoidal variations in airfoil surface.
+        Variation oscillates around the original shape.
+
+        s_start: float [m]
+            curve location where perturbation starts. Defined from LE.
+            +ve increasing towards suction side.
+        s_end: float [m]
+            curve location where perturbation ends. Defined from LE.
+            +ve increasing towards suction side.
+        act_len: float [m]
+            transition length between undisturbed and disturbed region.
+        amp: float [m]
+            amplitude of oscillation
+        lam: float [m]
+            wave length of oscillation
+        ni: integer [-] multiple of 8
+            additional grid points to resolve perturbed region
+        """
+        self.s_start = self.main_seg.sLE + s_start/self.main_seg.smax
+        self.s_end = self.main_seg.sLE + s_end/self.main_seg.smax
+        slen = self.s_end - self.s_start
+        # Change point distribution
+        main_ni = self.main_ni
+        main_seg = self.main_seg.copy()
+        i_start = np.argmin(abs(self.s_start-self.main_seg.s))
+        i_end = np.argmin(abs(self.s_end-self.main_seg.s))
+        i_LE = main_ni// 2 + 1
+        # Extra points are mostly concentrated in the oscillating area,
+        # but two eigth are placed outside the region as well
+        ni_o = ni//8
+        # Depending where the region is, some points need to be added
+        # at the LE sto leave the intial point distribution unaffected.
+        if self.s_start < self.main_seg.sLE:
+            i_LE += 4*ni_o
+        if self.s_end < self.main_seg.sLE:
+            i_LE += 4*ni_o
+        main_ni = main_ni + ni
+        i_start += ni_o
+        i_end += 7*ni_o
+        ds = slen/ni
+        # Increase resolution in perturbed area
+        dist = [[0., self.ds_TE/self.main_seg.smax, 1],
+                [self.s_start, ds/self.main_seg.smax, i_start+1],
+                [self.main_seg.sLE, self.ds_LE/self.main_seg.smax, i_LE+1],
+                [self.s_end, ds/self.main_seg.smax, i_end+1],
+                [1., self.ds_TE/self.main_seg.smax, main_ni]]
+        # Sort dist by the curve length
+        dist = sorted(dist, key=lambda x: x[0])
+        self.dist = dist
+        self.ds = ds*self.main_seg.smax
+
+        main_seg._compute_dp_flag = True
+        main_seg.redistribute(dist=dist)
+        # Scaling
+        lam /= main_seg.smax
+        act_len /= main_seg.smax
+        s = main_seg.s
+        # Perturbation in the curve normal direction
+        dn = amp/2*np.sin(2*np.pi*s/lam)
+        #print(dn)
+        # Activation function
+        tanc = 1./act_len
+        activ = (np.tanh(4*tanc*(s-self.s_start))+1)/2-(np.tanh(4*tanc*(s-self.s_end))+1)/2
+        # Final perturbation
+        dn *= activ
+        normals = np.array([-main_seg.dp[:,1],main_seg.dp[:,0]]).T
+
+        wavy = np.zeros(normals.shape[:])
+        for i in range(dn.shape[0]):
+            #print('normals[i,:]*dn[i]',normals[i,:]*dn[i])
+            wavy[i,:] = normals[i,:]*dn[i] + main_seg.points[i,:]
+        self.wavy = wavy
+        # Assemble profile
+        points = np.append(self.te_seg.points[self.te_isplit:], wavy[1:], axis=0)
+        points = np.append(points, self.te_seg.points[1:self.te_isplit+1], axis=0)
+        self.points = points
+        self.initialize(self.points)
+
+        return self
+
+    def step(self, s_start, s_end, h, ds_h, ni_step):
+        """
+        Add or remove material from the airfoil. This creates steps at the
+        end of the perturbed region. The difference between step and slot
+        lies in the positioning of the defect.  
+
+        s_start: float [m]
+            curve location where perturbation starts. Defined from LE.
+            +ve increasing towards suction side.
+        s_end: float [m]
+            curve location where perturbation ends. Defined from LE.
+            +ve increasing towards suction side.
+        h: float [m]
+            size of material addition (+ve) or removal (-ve)
+        ds_h: float [m]
+            grid size over h
+        ni_step: int [-]
+            additional grid points to resolve perturbed region
+        """
+        # Translate into normalised coordinates
+        self.s_step_start = self.main_seg.sLE+s_start/self.main_seg.smax
+        self.s_step_end = self.main_seg.sLE+s_end/self.main_seg.smax
+        # Find the number of grid points per step
+        self.ni_h = int(max(np.round(abs(h)/ds_h),np.round(abs(h)/self.ds_LE)))+1
+        self.ds_h = ds_h
+        # Additional points not used in the step are used over the airfoil
+        # (Make sure there are enough points to allow this)
+        ni_step = ni_step - 2*(self.ni_h-1)
+        # Find airfoil points closest to start and end
+        i_step_start = np.argmin(abs(self.s_step_start-self.main_seg.s))
+        i_step_end = np.argmin(abs(self.s_step_end-self.main_seg.s))
+        # Point ratio
+        ni_rat = np.array([i_step_start,
+                           i_step_end-i_step_start,
+                           self.main_ni-i_step_end])
+        ni_rat = ni_rat/self.main_ni
+        # Weight the point ratio (giving most weight to perturbed area)
+        ni_step_rat = np.array([1, 8, 1])*ni_rat
+        ni_step_rat /= sum(ni_step_rat)
+        ni_step_rat = np.round(ni_step*ni_step_rat)
+        # Ensure total no. of points does not exceed the no. specified
+        ni_step_rat[1] += (ni_step - sum(ni_step_rat))
+        # New no. of points
+        main_ni = self.main_ni + ni_step
+        main_seg = self.main_seg.copy()
+        # New point defintion
+        i_step_start = int(i_step_start+ni_step_rat[0])
+        i_step_end = int(i_step_end+sum(ni_step_rat[0:2]))
+
+        dist = [[0., self.ds_TE/self.main_seg.smax, 1],
+                [self.s_step_start, self.ds_h/self.main_seg.smax, i_step_start+1],
+                [self.main_seg.sLE, self.ds_LE/self.main_seg.smax, main_ni// 2 + 1],
+                [self.s_step_end, self.ds_h/self.main_seg.smax, i_step_end+1],
+                [1., self.ds_TE/self.main_seg.smax, main_ni]]
+        # Sort by curve length
+        dist = sorted(dist, key=lambda x: x[0])
+        main_seg.redistribute(dist=dist)
+        self.step = main_seg
+
+        # Apply disturbance normal to curve
+        main_seg._compute_dp()
+        normals = np.array([-main_seg.dp[:,1],main_seg.dp[:,0]]).T
+        step_points = main_seg.points[i_step_start:i_step_end+1,:]+ h*normals[i_step_start:i_step_end+1,:]
+        # Create the normal lines marching from original curve
+        step1 = np.zeros((self.ni_h,2))
+        step2 = np.zeros((self.ni_h,2))
+        for i in range(self.ni_h):
+            step1[i,:] = np.sign(h)*i*self.ds_h*normals[i_step_start,:] + main_seg.points[i_step_start,:]
+        for i in range(self.ni_h):
+            step2[i,:] = np.sign(h)*i*self.ds_h*normals[i_step_end,:] + main_seg.points[i_step_end,:]
+        # Assemble pertubed region
+        step_points = np.append(step1, step_points[1:], axis=0)
+        step_points = np.append(step_points[:-1], step2[::-1], axis=0)
+        step_points = np.append(main_seg.points[:i_step_start],step_points, axis=0)
+        step_points = np.append(step_points, main_seg.points[i_step_end+1:], axis=0)
+        # Assemble profile
+        points = np.append(self.te_seg.points[self.te_isplit:], step_points[1:], axis=0)
+        points = np.append(points, self.te_seg.points[1:self.te_isplit+1], axis=0)
+        self.i_step_start = i_step_start
+        self.i_step_end = i_step_end
+        self.points = points
+        self.initialize(self.points)
+
+        return self
+
+    def moveLE(self, x_start, dy):
+        """
+        Moves the LE up (dy > 0) or down (dy < 0) in the y direction
+
+        x_start: float [m]
+            x location where to start movement of LE
+        dy: float [m]
+            movement of LE in the y direction
+        """
+        points = self.points.copy()
+        x = (x_start - points[:,0])
+        # Range of x
+        dx = np.amax(x)
+        x /=dx
+        # Leave y value untouched downstream of x location
+        x[x<0] = 0
+        # cubic y perturbation growth towards LE
+        print(x.shape[:])
+        dypoints = x**2
+        dypoints*= dy
+        self.dypoints = dypoints
+        self.x = x
+
+        # Final perturbed airfoil
+        points[:,1] += dypoints
+        self.points = points
+        self.initialize(self.points)
+
+        return self
+
+    def flatsanding(self, s_mid, s_len, ni):
+        """
+        Introduces a flat region on airfoil
+
+        s_mid: float [m]
+            curve location in middle of perturbation. Defined from LE.
+            +ve increasing towards suction side.
+        s_len: float [m]
+            length of flat region
+        ni: int [-] (multiple of 4)
+            additional no. of points resolving feature
+        """
+        # Normalised start and end of flat region
+        self.s_start = self.main_seg.sLE+(s_mid-s_len/2)/self.main_seg.smax
+        self.s_end = self.main_seg.sLE+(s_mid+s_len/2)/self.main_seg.smax
+        main_ni = self.main_ni
+        main_seg = self.main_seg.copy()
+        # Point no. at feature ends
+        i_start = np.argmin(abs(self.s_start-self.main_seg.s))
+        i_end = np.argmin(abs(self.s_end-self.main_seg.s))
+        # Distribute additional points around the feature ends
+        i_LE = main_ni// 2 + 1
+        ni_q = ni//4
+        # Adjust the no. of points at the LE
+        if self.s_start < self.main_seg.sLE:
+            i_LE += 2*ni_q
+        if self.s_end < self.main_seg.sLE:
+            i_LE += 2*ni_q
+        # Add additional points to the original
+        main_ni = main_ni + ni
+        i_start += ni_q
+        i_end += 3*ni_q
+        # Grid size at the ends of the feature (either the same as at LE or evenly spaced)
+        ds = min(self.ds_LE,s_len/ni)
+        dist = [[0., self.ds_TE/self.main_seg.smax, 1],
+                [self.s_start, ds/self.main_seg.smax, i_start+1],
+                [self.main_seg.sLE, self.ds_LE/self.main_seg.smax, i_LE],
+                [self.s_end, ds/self.main_seg.smax, i_end+1],
+                [1., self.ds_TE/self.main_seg.smax, main_ni]]
+        # Sort by curve length
+        dist = sorted(dist, key=lambda x: x[0])
+        self.dist = dist
+        main_seg.redistribute(dist=dist)
+        # Flat region
+        flat = Line(main_seg.points[i_start,:], main_seg.points[i_end,:], i_end-i_start+1)
+        dist = [[0., ds/flat.smax, 1],
+                [1, ds/flat.smax, flat.ni]]
+        flat.redistribute(dist=dist)
+
+        # Assemble airfoil with TE closure
+        points = np.append(self.te_seg.points[self.te_isplit:], main_seg.points[1:i_start+1], axis=0)
+        points = np.append(points, flat.points[1:], axis=0)
+        points = np.append(points, main_seg.points[i_end+1:], axis=0)
+        points = np.append(points, self.te_seg.points[1:self.te_isplit+1], axis=0)
+
+        self.flat = flat
+        self.points = points
+        self.initialize(self.points)
+
+        return self
+
+    def rough_paras(self,s_len,h):
+        """
+        Computes the inputs to the roughness box for input.dat, for
+        symmetric roughness about the LE. It returns the edge points of the 
+        roughness box in a np.array self.box and the roughness length as
+        self.rough_len.
+
+        s_len: float [m]
+            length of roughness on each side of the LE
+        h: float [m]
+            roughness height of sandpaper (not Nikuradse)
+        """
+        self.s_start = self.main_seg.sLE - s_len/self.main_seg.smax
+        self.s_end = self.main_seg.sLE + s_len/self.main_seg.smax
+        i_start = np.argmin(abs(self.s_start-self.main_seg.s))
+        i_end = np.argmin(abs(self.s_end-self.main_seg.s))
+        points = self.main_seg.points.copy()
+        self.i = np.array([i_start, i_end])
+        xmax = (points[i_start,0]+points[i_end,0])/2.
+        c = points[:,0].max() - points[:,0].min()
+        xmin = points[:,0].min()
+        ymin = points[:,1].min()
+        ymax = points[:,1].max()
+        buf = 0.05*c
+        self.box = np.array([xmin-buf, xmax, ymin-buf, ymax+buf])
+        # The ratio between Nikuradse's roughness scale and actual roughness length for
+        # sandpaper like surfaces is 2 for a roughness concentration of about 30%
+        ks = 2.*h
+        # In Knopp's roughness model the roughness is given by
+        d0 = 0.03*ks
+        self.rough_len = d0
+
+        return self
+
+    def roughpatch_paras(self, s_step, s_len, h, enforce_dir=True):
+        """
+        Similar to rough_paras but now the start location of the roughness can be
+        set anywhere along the airfoil. Returns self.box_up and self.box_low 
+        instead of a single box for each side of the airfoil. If the roughness does
+        not cross the LE (always True if enforce_dir=True) then box_low=box_up.
+
+        s_step: float [m]
+            curve location where perturbation starts. Defined from LE.
+            +ve increasing towards suction side.
+        s_len: float [m]
+            length of roughness patch
+        h: float [m]
+            roughness height of sandpaper (not Nikuradse)
+        enforce_dir: bool 
+            If "True" the end point of the roughpatch is always placed towards the TE i.e.
+            roughness will not cross the LE. This is not the case when "False" then s_len
+            is just added to the start location and if s_step sufficiently close to the LE
+            and s_len long enough the LE is crossed by the patch. This is also why the output 
+            of this function is an upper and lower roughness box, one for the suction and the 
+            other for the pressure side. If the LE is not crossed both boxes are identical. 
+
+        """
+        # Intialize roughness box such that they do not influence the airfoil
+        self.box_up = np.array([5., 5.1, 5., 5.1])
+        self.box_low = np.array([5., 5.1, 5., 5.1])
+
+        # Translate into normalised coordinates
+        sstart = self.main_seg.sLE+s_step/self.main_seg.smax
+        if enforce_dir:
+            # place end point towards TE
+            send = self.main_seg.sLE+(s_step+np.sign(s_step)*s_len)/self.main_seg.smax
+        else:
+            send = self.main_seg.sLE+(s_step+s_len)/self.main_seg.smax
+        if sstart < send:
+            self.s_start = sstart
+            self.s_end = send
+        else:
+            self.s_start = send
+            self.s_end = sstart
+        # Point no. at feature ends
+        i_start = np.argmin(abs(self.s_start-self.main_seg.s))
+        i_end = np.argmin(abs(self.s_end-self.main_seg.s))
+        self.i = np.array([i_start, i_end])
+        i_LE = self.main_ni// 2 + 1
+        points = self.main_seg.points.copy()
+        c = points[:,0].max() - points[:,0].min()
+        xs = points[[i_start,i_end],0]
+        ys = points[[i_start,i_end],1]
+        # Determine the coordinates of the roughness boxes
+        # Patch on upper surface ony
+        if (self.s_start >= self.main_seg.sLE):
+            self.box_up = np.array([xs.min(), xs.max(), ys.min(), ys.max()])
+            self.box_low = self.box_up
+        elif (self.s_end <= self.main_seg.sLE):
+            self.box_low = np.array([xs.min(), xs.max(), ys.min(), ys.max()])
+            self.box_up = self.box_low
+        else:
+            self.box_up = np.array([points[:,0].min()-0.01*c, points[i_end,0],
+                                    points[i_LE,1], points[i_end,1]])
+
+            self.box_low = np.array([points[:,0].min()-0.01*c, points[i_start,0],
+                                    points[i_start,1], points[i_LE,1]])
+
+        # The ratio between Nikuradse's roughness scale and actual roughness length for
+        # sandpaper like surfaces is 2 for a roughness concentration of about 30%
+        ks = 2.*h
+        # In Knopp's roughness model the roughness is given by
+        d0 = 0.03*ks
+        self.rough_len = d0
+
+        return self
+
+    def smoothbite(self,bite_size):
+        """
+        Create over- (bite > 0) and underbite (bite < 0), which is then filled up 
+        again, modelling the effect of filling any imperfections at the LE. 
+
+        bite: float [m]
+            size of over or underbite
+        """
+        self.bite(bite_size)
+        i_LE = np.argmin(self.points[:,0])
+        i_ymax = np.argmax(self.points[:,1])
+        i_ymin = np.argmin(self.points[:,1])
+        pLE = self.points[i_LE]
+        # Vertical point abvoe/below LE point
+        pvLE = np.array([pLE[0], pLE[1] - 2*bite_size])
+
+        ang_p = np.arctan2(self.points[:,1]-pvLE[1],self.points[:,0]-pvLE[0])
+        if bite_size < 0:
+            ang_c = np.arctan2(self.dp[:,1],self.dp[:,0])
+            ang_diff = abs(ang_c[i_LE:i_ymax]-ang_p[i_LE:i_ymax])
+            i_tangent = np.argmin(ang_diff) + i_LE
+        else:
+            ang_c = np.arctan2(-self.dp[:,1],-self.dp[:,0])
+            ang_diff = abs(ang_c[i_ymin:i_LE]-ang_p[i_ymin:i_LE])
+            i_tangent = np.argmin(ang_diff) + i_ymin
+
+        self.i_LE = i_LE
+        self.i_tangent = i_tangent
+        self.pvLE = pvLE
+
+        # The filler surface is assumed to follow a Bezier Curve
+        c = BezierCurve()
+        c.ni = 100
+        c.add_control_point(np.array([pLE[0], pLE[1], 0]))
+        c.add_control_point(np.array([pvLE[0], pvLE[1], 0]))
+        c.add_control_point(np.array([self.points[i_tangent,0], self.points[i_tangent,1], 0]))
+        c.update()
+        self.c = c
+        # Depending on the bite sign,the curve needs to be build differently.
+        # The TE is excluded, as aftrwards the point distribution needs to be chnaged again
+        if bite_size < 0:
+            points = np.append(self.points[(self.ni_TE+1)//2:i_LE], c.points[1:,:2], axis=0)
+            points = np.append(points, self.points[i_tangent+1:-(self.ni_TE+1)//2], axis=0)
+        else:
+            points = np.append(self.points[(self.ni_TE+1)//2:i_tangent], c.points[-2::-1,:2], axis=0)
+            points = np.append(points, self.points[i_LE+1:-(self.ni_TE+1)//2], axis=0)
+
+        self.over_points = self.points
+        self.points = points
+        self.redistribute(self.ni-1,close_te=True)
+        self.initialize(self.points)
+
+        return self
+
+    def slot(self, s_step, s_len, h, ds_h, ni_step):
+        """
+        Add or remove material from the airfoil over a certain region.
+
+        s_step: float [m]
+            curve location where perturbation starts. Defined from LE.
+            +ve increasing towards suction side.
+        s_len: float [m]
+            length of perturbation 
+        h: float [m]
+            size of material addition (+ve) or removal (-ve)
+        ds_h: float [m]
+            grid size over h
+        ni_step: int [-]
+            additional grid points to resolve perturbed region
+        """
+        # Translate into normalised coordinates
+        sstart = self.main_seg.sLE+s_step/self.main_seg.smax
+        # Always place end point towards TE
+        send = self.main_seg.sLE+(s_step+np.sign(s_step)*s_len)/self.main_seg.smax
+        if sstart < send:
+            self.s_step_start = sstart
+            self.s_step_end = send
+        else:
+            self.s_step_start = send
+            self.s_step_end = sstart
+
+        # Find the number of grid points per step
+        self.ni_h = int(max(np.round(abs(h)/ds_h),np.round(abs(h)/self.ds_LE)))+1
+        self.ds_h = ds_h
+        # Additional points not used in the step are used over the airfoil
+        # (Make sure there are enough points to allow this)
+        ni_step = ni_step - 2*(self.ni_h-1)
+        # Find airfoil points closest to start and end
+        i_step_start = np.argmin(abs(self.s_step_start-self.main_seg.s))
+        i_step_end = np.argmin(abs(self.s_step_end-self.main_seg.s))
+        main_seg = self.main_seg.copy()
+        # Make sure all extra points are used on feature
+        i_step_start += ni_step//2
+        i_step_end += ni_step
+        i_LE = self.main_ni//2 + 1
+        if self.s_step_end < self.main_seg.sLE:
+            i_LE += ni_step
+
+        dist = [[0., self.ds_TE/self.main_seg.smax, 1],
+                [self.s_step_start, self.ds_h/self.main_seg.smax, i_step_start+1],
+                [self.main_seg.sLE, self.ds_LE/self.main_seg.smax, i_LE],
+                [self.s_step_end, self.ds_h/self.main_seg.smax, i_step_end+1],
+                [1., self.ds_TE/self.main_seg.smax, self.main_ni+ni_step]]
+        # Sort by curve length
+        dist = sorted(dist, key=lambda x: x[0])
+        main_seg.redistribute(dist=dist)
+        self.step = main_seg.points.copy()
+
+        # Apply disturbance normal to curve
+        main_seg._compute_dp()
+        normals = np.array([-main_seg.dp[:,1],main_seg.dp[:,0]]).T
+        step_points = main_seg.points[i_step_start:i_step_end+1,:]+ h*normals[i_step_start:i_step_end+1,:]
+        # Create the normal lines marching from original curve
+        step1 = np.zeros((self.ni_h,2))
+        step2 = np.zeros((self.ni_h,2))
+        for i in range(self.ni_h):
+            step1[i,:] = np.sign(h)*i*self.ds_h*normals[i_step_start,:] + main_seg.points[i_step_start,:]
+        for i in range(self.ni_h):
+            step2[i,:] = np.sign(h)*i*self.ds_h*normals[i_step_end,:] + main_seg.points[i_step_end,:]
+        # Assemble pertubed region
+        step_points = np.append(step1, step_points[1:], axis=0)
+        step_points = np.append(step_points[:-1], step2[::-1], axis=0)
+        step_points = np.append(main_seg.points[:i_step_start],step_points, axis=0)
+        step_points = np.append(step_points, main_seg.points[i_step_end+1:], axis=0)
+        # Assemble profile
+        points = np.append(self.te_seg.points[self.te_isplit:], step_points[1:], axis=0)
+        points = np.append(points, self.te_seg.points[1:self.te_isplit+1], axis=0)
+
+        self.i_step_start = i_step_start
+        self.i_step_end = i_step_end
+        self.points = points
+        self.initialize(self.points)
+
+        return self
+
+    def smoothslot_start(self, s_step, s_len, h, ds_h, ni_step):
+        """
+        Remove/add material with step at the start of the perturbation. 
+
+        s_step: float [m]
+            curve location where perturbation starts. Defined from LE.
+            +ve increasing towards suction side.
+        s_len: float [m]
+            length of perturbation 
+        h: float [m]
+            size of material addition (+ve) or removal (-ve)
+        ds_h: float [m]
+            grid size over h
+        ni_step: int [-]
+            additional grid points to resolve perturbed region
+        """
+        # A dditional class hosting all the LE modifications
+        # Translate into normalised coordinates
+        self.s_step_start = self.main_seg.sLE+s_step/self.main_seg.smax
+        # Always place end point towards TE
+        self.s_step_end = self.main_seg.sLE+(s_step+np.sign(s_step)*s_len)/self.main_seg.smax
+        # Find the number of grid points per step
+        self.ni_h = int(max(np.round(abs(h)/ds_h),np.round(abs(h)/self.ds_LE)))+1
+        self.ds_h = ds_h
+        # Additional points not used in the step are used over the airfoil
+        # (Make sure there are enough points to allow this)
+        ni_step = ni_step - (self.ni_h-1)
+        # Find airfoil points closest to start and end
+        i_step_start = np.argmin(abs(self.s_step_start-self.main_seg.s))
+        i_step_end = np.argmin(abs(self.s_step_end-self.main_seg.s))
+
+        if self.s_step_end < self.s_step_start:
+            i_step_end += ni_step//2
+            i_step_start += ni_step
+        else:
+            i_step_start += ni_step//2
+            i_step_end += ni_step
+        # New no. of points
+        main_ni = self.main_ni
+        main_seg = self.main_seg.copy()
+        # If the slot is before the LE than add some points to the LE
+        i_LE = main_ni//2 + 1
+        if self.s_step_start < self.main_seg.sLE:
+            i_LE += ni_step
+        # Coarsen towards the end of the step
+        ds_end = s_len/abs(h)*self.ds_h/self.main_seg.smax
+
+        dist = [[0., self.ds_TE/self.main_seg.smax, 1],
+                [self.s_step_start, self.ds_h/self.main_seg.smax, i_step_start+1],
+                [self.main_seg.sLE, self.ds_LE/self.main_seg.smax, i_LE],
+                [self.s_step_end, ds_end, i_step_end+1],
+                [1., self.ds_TE/self.main_seg.smax, main_ni+ni_step]]
+        # Sort by curve length
+        dist = sorted(dist, key=lambda x: x[0])
+        main_seg.redistribute(dist=dist)
+        self.step = main_seg.points.copy()
+
+        # Apply disturbance normal to curve
+        main_seg._compute_dp()
+        normals = np.array([-main_seg.dp[:,1],main_seg.dp[:,0]]).T
+        # Create the normal lines marching from original curve
+        step = np.zeros((self.ni_h,2))
+        for i in range(self.ni_h):
+            step[i,:] = np.sign(h)*i*self.ds_h*normals[i_step_start,:] + main_seg.points[i_step_start,:]
+        slen = s_len/self.main_seg.smax
+        if self.s_step_start < self.s_step_end:
+            step_reg = main_seg.points[i_step_start:i_step_end+1,:]
+            step_s = main_seg.s[i_step_start:i_step_end+1]
+            step_reg[:,0] += ((step_s-self.s_step_end)/slen)**2*h*normals[i_step_start:i_step_end+1,0]
+            step_reg[:,1] += ((step_s-self.s_step_end)/slen)**2*h*normals[i_step_start:i_step_end+1,1]
+            step_points = np.append(step,step_reg[1:,:], axis=0)
+            step_points = np.append(main_seg.points[:i_step_start],step_points, axis=0)
+            step_points = np.append(step_points, main_seg.points[i_step_end+1:], axis=0)
+        else:
+            step_reg = main_seg.points[i_step_end:i_step_start+1,:]
+            step_s = main_seg.s[i_step_end:i_step_start+1]
+            step_reg[:,0] += ((step_s-self.s_step_end)/slen)**2*h*normals[i_step_end:i_step_start+1,0]
+            step_reg[:,1] += ((step_s-self.s_step_end)/slen)**2*h*normals[i_step_end:i_step_start+1,1]
+            step_points = np.append(step_reg[:-1], step[::-1,:], axis=0)
+            step_points = np.append(main_seg.points[:i_step_end],step_points, axis=0)
+            step_points = np.append(step_points, main_seg.points[i_step_start+1:], axis=0)
+
+        # Assemble profile
+        points = np.append(self.te_seg.points[self.te_isplit:], step_points[1:], axis=0)
+        points = np.append(points, self.te_seg.points[1:self.te_isplit+1], axis=0)
+
+        self.i_step_start = i_step_start
+        self.i_step_end = i_step_end
+        self.points = points
+        self.initialize(self.points)
+
+        return self
+
+    def smoothslot_end(self, s_step, s_len, h, ds_h, ni_step):
+        """
+        Remove/add material with step at the end of the perturbation.  
+
+        s_step: float [m]
+            curve location where perturbation starts. Defined from LE.
+            +ve increasing towards suction side.
+        s_len: float [m]
+            length of perturbation 
+        h: float [m]
+            size of material addition (+ve) or removal (-ve)
+        ds_h: float [m]
+            grid size over h
+        ni_step: int [-]
+            additional grid points to resolve perturbed region
+        """
+        # Additional class hosting all the LE modifications
+        # Always place end point towards TE
+        self.s_step_end = self.main_seg.sLE+s_step/self.main_seg.smax
+        # Translate into normalised coordinates
+        self.s_step_start = self.main_seg.sLE+(s_step-np.sign(s_step)*s_len)/self.main_seg.smax
+        # Find the number of grid points per step
+        self.ni_h = int(max(np.round(abs(h)/ds_h),np.round(abs(h)/self.ds_LE)))+1
+        self.ds_h = ds_h
+        # Additional points not used in the step are used over the airfoil
+        # (Make sure there are enough points to allow this)
+        ni_step = ni_step - (self.ni_h-1)
+        # Find airfoil points closest to start and end
+        i_step_start = np.argmin(abs(self.s_step_start-self.main_seg.s))
+        i_step_end = np.argmin(abs(self.s_step_end-self.main_seg.s))
+        ds_orig = self.ds[i_step_start]
+        if self.s_step_end < self.s_step_start:
+            i_step_end += ni_step//2
+            i_step_start += ni_step
+        else:
+            i_step_start += ni_step//2
+            i_step_end += ni_step
+        # New no. of points
+        main_ni = self.main_ni
+        main_seg = self.main_seg.copy()
+        # If the slot is before the LE than add some points to the LE
+        i_LE = main_ni//2 + 1
+
+        if self.s_step_start < self.main_seg.sLE:
+            i_LE += ni_step//2
+        if self.s_step_end < self.main_seg.sLE:
+            i_LE += ni_step//2
+        # Coarsen towards the end of the step
+        ds_end = s_len/abs(h)*self.ds_h/self.main_seg.smax
+        ds_end = min(ds_orig,ds_end)
+        d = np.array([self.s_step_start-self.main_seg.sLE,self.s_step_end-self.main_seg.sLE])
+        if (np.sign(d[0]) == np.sign(d[1]))or(d[0]==0)or(d[1]==0):
+            dist = [[0., self.ds_TE/self.main_seg.smax, 1],
+                    [self.s_step_start, 0.66*self.ds_LE/self.main_seg.smax,  i_step_start+1],
+                    [self.s_step_end,  self.ds_h/self.main_seg.smax, i_step_end+1],
+                    [1., self.ds_TE/self.main_seg.smax, main_ni+ni_step]]
+        else:
+            dist = [[0., self.ds_TE/self.main_seg.smax, 1],
+                    [self.s_step_start, ds_end,  i_step_start+1],
+                    [self.main_seg.sLE, self.ds_LE/self.main_seg.smax, i_LE],
+                    [self.s_step_end,  self.ds_h/self.main_seg.smax, i_step_end+1],
+                    [1., self.ds_TE/self.main_seg.smax, main_ni+ni_step]]
+        # Sort by curve length
+        dist = sorted(dist, key=lambda x: x[0])
+        self.dist = dist
+        main_seg.redistribute(dist=dist)
+        self.step = main_seg.points.copy()
+
+        # Apply disturbance normal to curve
+        main_seg._compute_dp()
+        normals = np.array([-main_seg.dp[:,1],main_seg.dp[:,0]]).T
+        # Create the normal lines marching from original curve
+        step = np.zeros((self.ni_h,2))
+        for i in range(self.ni_h):
+            step[i,:] = np.sign(h)*i*self.ds_h*normals[i_step_end,:] + main_seg.points[i_step_end,:]
+        slen = s_len/self.main_seg.smax
+        if self.s_step_start < self.s_step_end:
+            step_reg = main_seg.points[i_step_start:i_step_end+1,:]
+            step_s = main_seg.s[i_step_start:i_step_end+1]
+            step_reg[:,0] += ((step_s-self.s_step_start)/slen)**2*h*normals[i_step_start:i_step_end+1,0]
+            step_reg[:,1] += ((step_s-self.s_step_start)/slen)**2*h*normals[i_step_start:i_step_end+1,1]
+            step_points = np.append(step_reg[:-1,:],step[::-1,:], axis=0)
+            step_points = np.append(main_seg.points[:i_step_start],step_points, axis=0)
+            step_points = np.append(step_points, main_seg.points[i_step_end+1:], axis=0)
+        else:
+            step_reg = main_seg.points[i_step_end:i_step_start+1,:]
+            step_s = main_seg.s[i_step_end:i_step_start+1]
+            step_reg[:,0] += ((step_s-self.s_step_start)/slen)**2*h*normals[i_step_end:i_step_start+1,0]
+            step_reg[:,1] += ((step_s-self.s_step_start)/slen)**2*h*normals[i_step_end:i_step_start+1,1]
+            step_points = np.append(step[:-1,:],step_reg, axis=0)
+            step_points = np.append(main_seg.points[:i_step_end],step_points, axis=0)
+            step_points = np.append(step_points, main_seg.points[i_step_start+1:], axis=0)
+
+        # Assemble profile
+        points = np.append(self.te_seg.points[self.te_isplit:], step_points[1:], axis=0)
+        points = np.append(points, self.te_seg.points[1:self.te_isplit+1], axis=0)
+
+        self.i_step_start = i_step_start
+        self.i_step_end = i_step_end
+        self.points = points
+        self.initialize(self.points)
+
+        return self
+
+    def stallstrip(self, s_strip, h, ni, ds, ang):
+        """
+        Creates a triangular stallstrip with max height at s_strip.
+
+        s_strip: float [m]
+            curve location of strip. Defined from LE.
+            +ve increasing towards suction side.
+        h: float [m]
+            height of stall strip
+        ni: int [-] (multiple of 4)
+            additional no. of points resolving feature
+        ds: float [m]
+            grid size over step
+        ang: float [deg]
+            angle in tip corner
+        """
+        # The stall strip has triangular shape with equal side lenghts
+        L = h/np.cos(np.deg2rad(ang/2))
+        L_base = 2*h*np.tan(np.deg2rad(ang/2))
+        # Normalised start and end of flat region
+        self.s_start = self.main_seg.sLE + (s_strip-L_base/2)/self.main_seg.smax
+        self.s_mid = self.main_seg.sLE + s_strip/self.main_seg.smax
+        self.s_end = self.main_seg.sLE + (s_strip+L_base/2)/self.main_seg.smax
+        main_ni = self.main_ni
+        main_seg = self.main_seg.copy()
+        # Point no. at feature ends
+        i_start = np.argmin(abs(self.s_start-self.main_seg.s))
+        i_mid = np.argmin(abs(self.s_mid-self.main_seg.s))
+        i_end = np.argmin(abs(self.s_end-self.main_seg.s))
+        # Distribute additional points around the feature ends
+        i_LE = main_ni// 2 + 1
+        # Grid step on strip
+        ds_tr = min(self.ds_LE,ds)
+        # Number of points per side of triangle
+        ni_L = 2*(np.round(2*L/ds_tr)//4)
+        # Ensure not too many points will be used
+        if ni_L > 2*(ni//6):
+            ni_L = 2*(ni//6)
+            ds_tr = min(L/ni_L,self.ds_LE)
+        # Determine no. of points used around stall strip
+        ni_extra = (ni - 2.*ni_L)//2
+        ni_L = int(ni_L)
+        ni_extra = int(ni_extra)
+        ni_extra += (ni - 2*ni_extra - 2*ni_L)
+        # Adjust the no. of points at the LE
+        if self.s_start < self.main_seg.sLE:
+            i_add = int(ni_extra*min((self.main_seg.sLE-self.s_start)*self.main_seg.smax/(4*L),1.))
+            i_LE +=   i_add + ni_L//2
+        if self.s_mid < self.main_seg.sLE:
+            i_LE += ni_L
+        if self.s_end < self.main_seg.sLE:
+            i_add = int(ni_extra*min((self.main_seg.sLE-self.s_start)*self.main_seg.smax/(4*L),1.))
+            i_LE += (ni_L//2 + i_add)
+        # Add additional points to the original
+        main_ni = main_ni + ni
+        i_start += ni_extra
+        i_mid = (i_start+ ni_L)
+        i_end = (i_mid + ni_L)
+
+        if np.sqrt((self.s_mid-self.main_seg.sLE)**2)*self.main_seg.smax < (10*L):
+            dist = [[0., self.ds_TE/self.main_seg.smax, 1],
+                    [self.s_start, ds_tr/self.main_seg.smax, i_start+1],
+                    [self.s_mid, ds_tr/self.main_seg.smax, i_mid+1],
+                    [self.s_end, ds_tr/self.main_seg.smax, i_end+1],
+                    [1., self.ds_TE/self.main_seg.smax, main_ni]]
+        else:
+            dist = [[0., self.ds_TE/self.main_seg.smax, 1],
+                    [self.s_start, ds_tr/self.main_seg.smax, i_start+1],
+                    [self.s_mid, ds_tr/self.main_seg.smax, i_mid+1],
+                    [self.main_seg.sLE, self.ds_LE/self.main_seg.smax, i_LE],
+                    [self.s_end, ds_tr/self.main_seg.smax, i_end+1],
+                    [1., self.ds_TE/self.main_seg.smax, main_ni]]
+        # Sort by curve length
+        dist = sorted(dist, key=lambda x: x[0])
+        main_seg.redistribute(dist=dist)
+
+        # Find the position of the tip of the trip strip
+        main_seg._compute_dp()
+        normal = np.array([-main_seg.dp[i_mid,1],main_seg.dp[i_mid,0]]).T
+        tip_xy = normal*h+main_seg.points[i_mid,:]
+
+        l1 = Line(main_seg.points[i_start,:], tip_xy, i_mid-i_start+1)
+        l2 = Line(tip_xy, main_seg.points[i_end,:], i_end-i_mid+1)
+        l1.redistribute(dist= [[0., ds_tr/l1.smax, 1],
+                               [1., ds_tr/l1.smax, i_mid-i_start+1]])
+        l2.redistribute(dist=[[0., ds_tr/l2.smax, 1],
+                              [1., ds_tr/l2.smax, i_end-i_mid+1]])
+        tr = np.append(l1.points[1:], l2.points[1:], axis=0)
+
+        # Assemble airfoil with TE closure
+        points = np.append(self.te_seg.points[self.te_isplit:], main_seg.points[1:i_start+1], axis=0)
+        points = np.append(points, tr, axis=0)
+        points = np.append(points, main_seg.points[i_end+1:], axis=0)
+        points = np.append(points, self.te_seg.points[1:self.te_isplit+1], axis=0)
+
+        self.tr = tr
+        self.points = points
+        self.initialize(self.points)
+
+        return self            
 
 class BlendAirfoilShapes(object):
     """
diff --git a/PGL/test/test_airfoil.py b/PGL/test/test_airfoil.py
index 1b4cf4761d22e30dfd4d63aea2a46239a7dc56bb..3aedcf829bc0f663c09021cfa3eb5ca598567e2f 100644
--- a/PGL/test/test_airfoil.py
+++ b/PGL/test/test_airfoil.py
@@ -120,6 +120,105 @@ class AirfoilShapeTests2(unittest.TestCase):
         s = np.sum(af.points[:,1])
         self.assertAlmostEqual(s,0.6714187425159088, places=5)
 
+    def test_wavy(self):
+
+        af = AirfoilShape(points=np.loadtxt('data/ffaw3241.dat'))
+        af.redistribute(128, close_te=True)
+        af.wavy(-100e-3, 100e-3, 40e-3, 0.5e-3, 20e-3, 128)
+        s = np.sum(af.points[:,1])
+        self.assertAlmostEqual(s,1.0209455387348676, places=5)
+
+    def test_bite(self):
+
+        af = AirfoilShape(points=np.loadtxt('data/ffaw3241.dat'))
+        af.redistribute(128, close_te=True)
+        af.bite(-20e-3)
+        s = np.sum(af.points[:,1])
+        self.assertAlmostEqual(s,0.7387984401323768, places=5)
+        af.bite(20e-3)
+        s = np.sum(af.points[:,1])
+        self.assertAlmostEqual(s,0.6685177054337146, places=5)
+
+    def test_step(self):
+
+        af = AirfoilShape(points=np.loadtxt('data/ffaw3241.dat'))
+        af.redistribute(128, close_te=True)
+        af.step(-100e-3, 100e-3, 1e-3, 0.1e-3/3, 128)
+        s = np.sum(af.points[:,1])
+        self.assertAlmostEqual(s,0.9917633146652349, places=5)
+
+    def test_moveLE(self):
+
+        af = AirfoilShape(points=np.loadtxt('data/ffaw3241.dat'))
+        af.redistribute(128, close_te=3)
+        af.moveLE(0.15, 15e-3)
+        s = np.sum(af.points[:,1])
+        self.assertAlmostEqual(s,1.4250613036049513, places=5)
+
+    def test_flatsanding(self):
+
+        af = AirfoilShape(points=np.loadtxt('data/ffaw3241.dat'))
+        af.redistribute(128, close_te=True)
+        af.flatsanding(-10e-3,  25e-3, 128)
+        s = np.sum(af.points[:,1])
+        self.assertAlmostEqual(s,0.18489494967820264, places=5)
+
+    def test_rough_paras(self):
+
+        af = AirfoilShape(points=np.loadtxt('data/ffaw3241.dat'))
+        af.redistribute(128, close_te=3)
+        af.rough_paras(100e-3, 1e-3)
+        box_ref = [-0.05001856,  0.06530538, -0.16562343,  0.17532873]
+        self.assertEqual(np.testing.assert_allclose(af.box, box_ref, 1E-06), None)
+
+    def test_roughpatch_paras(self):
+
+        af = AirfoilShape(points=np.loadtxt('data/ffaw3241.dat'))
+        af.redistribute(128, close_te=3)
+        af.roughpatch_paras(10e-3, 100e-3, 1e-3)
+        box_ref = [0.00099846, 0.06610987, 0.00825168, 0.07490293]
+        self.assertEqual(np.testing.assert_allclose(af.box_up, box_ref, 1E-06), None)
+
+    def test_smoothbite(self):
+
+        af = AirfoilShape(points=np.loadtxt('data/ffaw3241.dat'))
+        af.redistribute(128, close_te=3)
+        af.smoothbite(20e-3)
+        s = np.sum(af.points[:,1])
+        self.assertAlmostEqual(s, 0.8913207863514554, places=5)
+
+    def test_slot(self):
+
+        af = AirfoilShape(points=np.loadtxt('data/ffaw3241.dat'))
+        af.redistribute(128, close_te=3)
+        af.slot(-20e-3, 5e-3, 2e-3, 0.1e-3/3, 128)
+        s = np.sum(af.points[:,1])
+        self.assertAlmostEqual(s,-1.6577561460630772, places=5)
+
+    def test_smoothslot_start(self):
+
+        af = AirfoilShape(points=np.loadtxt('data/ffaw3241.dat'))
+        af.redistribute(128, close_te=3)
+        af.smoothslot_start(-20e-3, 20e-3, 2e-3, 0.1e-3/3, 128)
+        s = np.sum(af.points[:,1])
+        self.assertAlmostEqual(s,-2.8872653934553076, places=5)
+
+    def test_smoothslot_end(self):
+
+        af = AirfoilShape(points=np.loadtxt('data/ffaw3241.dat'))
+        af.redistribute(128, close_te=3)
+        af.smoothslot_end(-20e-3, 20e-3, 2e-3, 0.1e-3/3, 128)
+        s = np.sum(af.points[:,1])
+        self.assertAlmostEqual(s,-1.8455353600696762, places=5)
+
+    def test_stallstrip(self):
+
+        af = AirfoilShape(points=np.loadtxt('data/ffaw3241.dat'))
+        af.redistribute(128, close_te=3)
+        af.stallstrip(0., 3e-3, 128, 0.1e-3/5, 90.)
+        s = np.sum(af.points[:,1])
+        self.assertAlmostEqual(s,1.1306133840208155, places=5)
+
 if __name__ == '__main__':
 
     unittest.main()