#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Package SDPMgeometry contains functions to create grids for the Source and 
Doublet Panel Method.

This code is part of the SDPMflut Matlab distribution.
Copyright (C) 2024 Grigorios Dimitriadis

This program is free software: you can redistribute it and/or modify it 
under the terms of the GNU General Public License as published by the 
Free Software Foundation, either version 3 of the License, or (at your 
option) any later version.
 
This program is distributed in the hope that it will be useful, but 
WITHOUT ANY WARRANTY; without even the implied warranty of 
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 
Public License for more details.

You should have received a copy of the GNU General Public License along 
with this program. If not, see <https://www.gnu.org/licenses/>.
"""

import ctypes
import numpy as np
import numpy.matlib
from numpy.ctypeslib import ndpointer
import platform
from scipy import linalg
from scipy.interpolate import CubicSpline
import airfoils
        
def normalarea(x,y,z):
    # Calculates the unit normal vector and surface area of a near-planar 
    # quadrilateral panel. The vertices of the quadrilateral are stored in 
    # vectors x, y and z. Also calculates a measure of the coplanarity of the 
    # four vertices.
    # x: 4-element array of the x-coordinates of the vertices of the panel
    # y: 4-element array of the y-coordinates of the vertices of the panel
    # z: 4-element array of the z-coordinates of the vertices of the panel
    # normal: 3-lement array of the unit vector normal to the panel
    # s: area of the panel
    # cpln: Measure of the coplanarity of the four vertices.

    # Normal vector from diagonals 13 and 24. The vector is normal to these two
    # diagonals. Equation 5.90
    normal=np.cross(np.array([x[2]-x[0], y[2]-y[0], z[2]-z[0]]),np.array([x[1]-x[3], y[1]-y[3], z[1]-z[3]]))
    normal=normal/linalg.norm(normal)
    # Calculate a second normal from edges 12 and 13. The vector is normal to 
    # these two edges.
    normal2=np.cross(np.array([x[2]-x[0], y[2]-y[0], z[2]-z[0]]),np.array([x[1]-x[0], y[1]-y[0], z[1]-z[0]]))
    mag2=linalg.norm(normal2)
    # Test if the normal and normal2 are colinear. If they are, their cross 
    # product will be equal to zero. Store the magnitude of the cross product 
    # in cpln.
    if mag2 != 0.0:
        normal2=normal2/mag2
        cpln=np.cross(normal,normal2)
        cpln=linalg.norm(cpln)
    else:
        cpln=0.0
    # Enf if
    # Calculate a third normal from diagonal 13 and edge 14. The vector is 
    # normal to these two segments.
    normal3=np.cross(np.array([x[2]-x[0], y[2]-y[0], z[2]-z[0]]),np.array([x[3]-x[0], y[3]-y[0], z[3]-z[0]]))
    mag3=linalg.norm(normal3)
    # Calculate panel area from equation 5.91
    s=mag2/2.0+mag3/2.0;
    
    return normal,s,cpln

def SDPMgeometry_trap_fun(body,ibody,mhalf,nhalf,nchords,Dxw,mirroredwing,linchord,linspan,trap,name,dir_tau,rollpitchyaw,rollpitchyaw_cent,lexyz,nmin,cref):
    # Calculates the coordinates of the panel vertices for the Source and
    # Doublet Panel method.    
    # body: struct array containing the geometry of all the bodies in the flow.
    #       Each element of body describes a different wing, fairing or 
    #       fuselage
    # ibody: Index of body for which to determine the self-influence
    # mhalf: number of chordwise wing panels on one side of wing (upper or
    #         lower).
    # nchords: Length of wake in mean chords.
    # Dxw: Chordwise wake spacing parameter. Dxw=0: No wake, Dxw=-1: One 
    #        chordwise panel, Dxw>0: The chordwise wake spacing will be 
    #        Dxw times the mean chordwise spacing of the corresponding 
    #        wing.
    # nhalf: number of spanwise panels on half-wing. This may change if there
    #        are more than one trapezoidal section and some of them have less  
    #        than nmin panels    
    # mirroredwing: If mirroredwing=-1: a left half-wing will be created
    #               If mirroredwing=1: a right half-wing will be created
    #               If mirroredwing=2: two mirrored half-wings will be created.
    # linchord: Chordwise panel distribution: linchord=1 constant, 
    #           linchord=2 denser at the leading edge. 
    # linspan: Spanwise panel distribution: linspan=1 constant, linspan=2 
    #           denser at the wing tip(s) and midspan. 
    # trap: struct array containing wing geometry information arranged in 
    #           trapezoidal sections.
    # name: String of the name of the wing.
    # dir_tau: direction in which the unit tangent vector  for this wing
    #           (tauxx, tauxy, tauxz) has a zero component. Can be 1, 2 or 3.
    # rollpitchyaw: 3-element array containing rotation angles around the x, y
    #           y and z axes.
    # rollpitchyaw_cent: 3-element of the x, y and z coordinates around which
    #           the rotations in rollpitchyaw are to be applied.
    # lexyz: 3-element of the x, y and z coordinates of the root leading edge
    #           of the wing.
    # nmin: Minimum number of spanwise panels per trapezoidal section.
    # cref: Reference chord length (usually mean aerodynamic chord)
    
    # Total number of chordwise panels (upper and lower sides)
    m=2*mhalf

    # Calculate half span
    bhalf=np.sum(trap['span'])
    # Incremental array of half spans
    bsum=np.zeros(len(trap) + 1)
    bsum[1:len(trap) + 1]=np.cumsum(trap['span'])

    # Calculate the planform area of each trapezoidal section
    Strap=(trap['rootchord']+trap['rootchord']*trap['taper'])*trap['span']/2.0
    # Calculate wing planform area
    S=np.sum(Strap)
    # Calculate wing aspect ratio
    AR=(2.0*bhalf)**2/(2.0*S)

    # Define chordwise panel grid 
    if linchord == 1:
        # Linearly spaced non-dimensional chordwise panel vertex coordinates
        xp=np.linspace(0,1,num=mhalf+1) 
        xp=np.concatenate([np.flip(xp),xp[1:m+1]])
    else:
        # Nonlinearly spaced non-dimensional chordwise panel vertex coordinates from equation 5.182
        xp=1-np.sin(np.linspace(0,np.pi,num=m+1))
    # End if
    
    # Distribute spanwise panels over the trapezoidal sections
    n_sec=np.round(trap['span']/bhalf*nhalf)
    n_sec=n_sec.astype(int)
    # Ensure that none of the sections have fewer than nmin spanwise panels
    iko=np.argwhere(n_sec <= nmin)
    if len(iko) > 0:
        n_sec[iko]=nmin
    # End if
    # Update spanwise panel number
    nhalf=np.sum(n_sec)
 
    if mirroredwing == 2:
        n=2*nhalf # Number of spanwise panels for the full wing
    else:
        n=nhalf
    # end if
    
    # Initialize grid
    Xp0=np.zeros((m+1,nhalf+1))
    Yp0=np.zeros((m+1,nhalf+1))
    Zp0=np.zeros((m+1,nhalf+1))
    
    # Cumulative sum of panels
    nsum=np.zeros(len(trap) + 1,dtype=int)
    nsum[1:len(trap) + 1]=np.cumsum(n_sec)
    # Set root x-coordinate of first section to zero
    xrootle=0.0
    # Cycle over the trapezoidal sections
    for i in range(0,len(trap)): 
        # Set up spanwise panel grid
        if n_sec[i] == 0:
            # numpy.linspace outputs the start value when num=1, whereas Matlab's linspace outputs the stop value.
            yp=np.array([1.0])
        else:
            if linspan == 1:
                # Linearly spaced non-dimensional spanwise panel vertex coordinates
                yp=np.linspace(0,1,num=n_sec[i]+1) 
            else:
                # Nonlinearly spaced non-dimensional spanwise panel vertex coordinates from equation 5.183
                theta=np.linspace(np.pi,0,num=n_sec[i]+1)
                yp=(1+np.cos(theta))/2.0
            # End if
        # End if        
        # Calculate linear variation of chord length over the trapezoidal section
        cvar=(trap['taper'][i]-1)*trap['rootchord'][i]*np.absolute(yp)+trap['rootchord'][i]
        # Calculate leading edge x-coordinate variation along the section due
        # to sweep
        xsweep=np.tan(trap['sweepLE'][i])*abs(yp)*trap['span'][i]
        # Calculate x-coordinates of panel vertices
        Xp0[:,nsum[i]:nsum[i+1]+1]=np.transpose(np.matlib.repmat(xp,n_sec[i]+1,1))*np.matlib.repmat(cvar,m+1,1)+\
            np.matlib.repmat(xsweep+xrootle,m+1,1)+trap['xledist'][i]
        # Calculate y-coordinates of panel vertices
        Yp0[:,nsum[i]:nsum[i+1]+1]=np.matlib.repmat(yp,m+1,1)*trap['span'][i]+bsum[i]
 
        # Calculate non-dimensional root airfoil shape
        zp1,zpcamb=airfoils.get_airfoil_shape(trap['rootairfoil'][i],xp,mhalf,trap['rootairfoilparams'][i])        
        # Calculate non-dimensional tip airfoil shape  
        zp2,zpcamb=airfoils.get_airfoil_shape(trap['tipairfoil'][i],xp,mhalf,trap['tipairfoilparams'][i])
     
        # Interpolate linearly between the root and tip airfoils  
        zpmat=np.zeros((m+1,n_sec[i]+1))
        for iy in range(0,n_sec[i]+1): 
            zpmat[:,iy]=(zp2-zp1)*yp[iy]+zp1
        # End for
        
        # Calculate z-coordinates of panel vertices
        if trap['rootairfoil'][i] == 'flatplate':
            Zp0[:,nsum[i]:nsum[i+1]+1]=zpmat*trap['tipairfoilparams'][i][0]
        else:
            Zp0[:,nsum[i]:nsum[i+1]+1]=zpmat*np.matlib.repmat(cvar,m+1,1)
        # End for
        # Calculate  root x-coordinate of next section
        xrootle=xrootle+xsweep[n_sec[i]]+trap['xledist'][i]
    # End for
    
    # Only executed for flat plate airfoils
    if trap['rootairfoil'][0] == 'flatplate':
        # We assume that if the first trapezoidal section is a flat plate, the
        # entire wing is a flat plate of the same thickness. We then round off the
        # leading edge and bevel the trailing edge.
        # Bevel angle of airfoil (rounded to nearest degree)
        bangle=round(np.arctan(trap['rootairfoilparams'][0][0]/trap['rootairfoilparams'][0][1])*180.0/np.pi)*np.pi/180.0
        # Round off leading edge
        for i in range(0,nhalf+1):
            # Determine x coordinates of points that lie in the rounded off
            # section    
            iko=np.argwhere(Xp0[:,i]-Xp0[mhalf,i] <= trap['rootairfoilparams'][0][0])
            # Radius of circle
            rcirc=trap['rootairfoilparams'][0][0]
            # Centre of circle
            xcent=Xp0[mhalf,i]+rcirc
            zcent=0.0
            # Circle equation
            zz=zcent+np.sqrt(rcirc**2-(Xp0[iko,i]-xcent)**2+1e-16) # Add 1e-16 to avoid tiny negative argument in square root
            # Calculate z coordinates of points lying on a circle centred on
            # the chordline, one airfoil thickness behind the leading edge.
            Zp0[iko,i]=zz*np.sign(Zp0[iko,i])
        # End for
        # Bevel the trailing edge
        for i in range(0,nhalf+1):
            # Determine x coordinates of points that lie in the beveled
            # section            
            iko=np.argwhere(Xp0[0,i]-Xp0[:,i] <= trap['rootairfoilparams'][0][1])
            if len(iko) == 2:
                # There are no such points except for the trailing edge.
                # Move the x coordinates of the penultimate trailing edge
                # points to the edge of the bevel.
                Xp0[1,i]=Xp0[0,i]-trap['rootairfoilparams'][0][1]
                Xp0[m-1,i]=Xp0[0,i]-trap['rootairfoilparams'][0][1]
            else:
                # Calculate z coordinates of points lying on a straight line
                # passing from the trailing edge and parallel to the bevel
                # angle. 
                dummy=-(Xp0[iko,i]-Xp0[0,i])*np.tan(bangle)*np.sign(iko-mhalf)+Zp0[0,i]
                # Only apply bevel at a point whose current Z coordinate is higher 
                # than the bevelled Z coordinate.
                dd=np.concatenate([Zp0[iko,i],dummy],axis=1)
                Zp0[iko,i]=np.expand_dims(np.min(np.abs(dd),axis=1),axis=1)*np.sign(Zp0[iko,i])
            # End if
        # End for
    # End if
        
    # Apply the twist angle variation
    for i in range(0,len(trap)): 
        # Set up spanwise panel grid
        if linspan == 1:
            # Linearly spaced non-dimensional spanwise panel vertex coordinates
            yp=np.linspace(0,1,num=n_sec[i]+1) 
        else:
            # Nonlinearly spaced non-dimensional spanwise panel vertex coordinates from equation 5.183
            theta=np.linspace(np.pi,0,num=n_sec[i]+1)
            yp=(1+np.cos(theta))/2.0
        # End if        
        # Calculate linear variation of chord length over the trapezoidal section
        cvar=(trap['taper'][i]-1)*trap['rootchord'][i]*np.absolute(yp)+trap['rootchord'][i]
        # Calculate linear twist variation from root to tip of trapezoidal
        # section
        epsilonvar=trap['roottwist'][i]+(trap['tiptwist'][i]-trap['roottwist'][i])*np.absolute(yp)
        # For the first trapezoidal section, rotate all spanwise sections. For
        # all other trapezoidla sections, do not rotate the first spanwise section
        # because it has already been rotated
        if i == 1:
            istart=0
        else:
            istart=1
        # End if
        # Rotate the airfoils
        for iy in range(istart,n_sec[i]+1): 
            # Calculate rotation matrix for the current spanwise station
            Ry=np.array([[np.cos(epsilonvar[iy]), np.sin(epsilonvar[iy])],[-np.sin(epsilonvar[iy]), np.cos(epsilonvar[iy])]])
            # Calculate rotation centre
            xyzcent=np.array([Xp0[mhalf,nsum[i]+iy], Yp0[mhalf,nsum[i]+iy], Zp0[mhalf,nsum[i]+iy]])+\
                trap['twistcent'][i]*np.array([cvar[iy], 0.0, 0.0])
            # Apply twist rotation around twistcent at the current spanwise station
            dummy=Ry @ np.array([Xp0[:,nsum[i]+iy]-xyzcent[0],Zp0[:,nsum[i]+iy]-xyzcent[2]])
            # Assign rotated coordinates
            Xp0[:,nsum[i]+iy]=dummy[0,:]+xyzcent[0]    # x-coordinates of panel vertices
            Zp0[:,nsum[i]+iy]=dummy[1,:]+xyzcent[2]    # z-coordinates of panel vertices
        # End for
    # End for

    # Calculate non-dimensional spanwise coordinates of all wing panel vertices
    yp=Yp0[0,:]/bhalf

    # Dihedral calculation
    # Calculate the difference in dihedral angle between sections
    if len(trap) == 1:
        dihedral_diff=np.array(trap['dihedral'][0],ndmin=1)
    else:
        dihedral_diff=np.concatenate([np.array([trap['dihedral'][0]]), np.diff(trap['dihedral'])])  
    # End if
    # Apply dihedral angle
    for i in range(0,len(trap)): 
        # Centre of rotation (root leading edge of trapezoidal section)
        xyzcent=np.array([Xp0[mhalf,nsum[i]], Yp0[mhalf,nsum[i]], Zp0[mhalf,nsum[i]]])
        # Rotation matrix
        Rx=np.array([[1, 0, 0],[0, np.cos(dihedral_diff[i]), -np.sin(dihedral_diff[i])], \
                     [0, np.sin(dihedral_diff[i]), np.cos(dihedral_diff[i])]])
        # Apply dihedral rotation around root leading edge of trapezoidal section
        xx=np.reshape(Xp0[:,nsum[i]+1:nhalf+1],(1,(m+1)*(nhalf-nsum[i])),order='C')-xyzcent[0]
        yy=np.reshape(Yp0[:,nsum[i]+1:nhalf+1],(1,(m+1)*(nhalf-nsum[i])),order='C')-xyzcent[1]
        zz=np.reshape(Zp0[:,nsum[i]+1:nhalf+1],(1,(m+1)*(nhalf-nsum[i])),order='C')-xyzcent[2]
        dummy=Rx @ np.concatenate((xx,yy,zz),axis=0)
        # Assign rotated coordinates
        Xp0[:,nsum[i]+1:nhalf+1]=np.reshape(dummy[0,:],(m+1,nhalf-nsum[i]),order='C')+xyzcent[0]
        Yp0[:,nsum[i]+1:nhalf+1]=np.reshape(dummy[1,:],(m+1,nhalf-nsum[i]),order='C')+xyzcent[1]
        Zp0[:,nsum[i]+1:nhalf+1]=np.reshape(dummy[2,:],(m+1,nhalf-nsum[i]),order='C')+xyzcent[2]
        if i > 0:
            # Also rotate the root of the current trapezoidal section (tip of
            # previous trapezoidal section) by half the dihedral angle
            # difference.
            # Rotation matrix
            Rx=np.array([[1, 0, 0],[0, np.cos(dihedral_diff[i]/2.0), -np.sin(dihedral_diff[i]/2.0)], \
                        [0, np.sin(dihedral_diff[i]/2.0), np.cos(dihedral_diff[i]/2.0)]])
            # Apply dihedral rotation around root leading edge of trapezoidal section
            xyz=np.zeros((3,m+1))
            xyz[0,:]=Xp0[:,nsum[i]]-xyzcent[0]
            xyz[1,:]=Yp0[:,nsum[i]]-xyzcent[1]
            xyz[2,:]=Zp0[:,nsum[i]]-xyzcent[2]
            dummy=Rx @ xyz
            # Assign rotated coordinates
            Xp0[:,nsum[i]]=dummy[0,:]+xyzcent[0]
            Yp0[:,nsum[i]]=dummy[1,:]+xyzcent[1]
            Zp0[:,nsum[i]]=dummy[2,:]+xyzcent[2]
    # End for
    
    # Re-arrange wing panels if necessary 
    if mirroredwing == -1:
        # Left half-wing has been requested, mirror panel vertices
        Xp0=np.flip(Xp0,axis=1)
        Yp0=-np.flip(Yp0,axis=1)
        yp=-np.flip(yp)
        Zp0=np.flip(Zp0,axis=1)
    elif mirroredwing == 2:
        # Full wing has been requested, mirror panel vertices and concatenate with original
        if nhalf == 0:
            Xp0=np.concatenate([Xp0,Xp0],axis=1)
            Yp0=np.concatenate([-Yp0,Yp0],axis=1)
            yp=np.concatenate([-yp,yp],axis=0)
            Zp0=np.concatenate([Zp0,Zp0],axis=1)
            n=1
        else:
            Xp0=np.concatenate([np.flip(Xp0,axis=1),Xp0[:,1:nhalf+1]],axis=1)
            Yp0=np.concatenate([-np.flip(Yp0,axis=1),Yp0[:,1:nhalf+1]],axis=1)
            yp=np.concatenate([-np.flip(yp),yp[1:nhalf+1]],axis=0)
            Zp0=np.concatenate([np.flip(Zp0,axis=1),Zp0[:,1:nhalf+1]],axis=1)
        # End if
    # End if
    
    # Move entire wing to required leading edge position
    Xp0=Xp0+lexyz[0]
    Yp0=Yp0+lexyz[1]
    Zp0=Zp0+lexyz[2]

    # Roll the entire wing
    if rollpitchyaw[0] != 0.0 :
        # Rotation angle
        rollangle=rollpitchyaw[0]
        # Rotation centre
        xf=rollpitchyaw_cent[0]
        yf=rollpitchyaw_cent[1]
        zf=rollpitchyaw_cent[2]
        # Rotation matrix
        Rx=np.array([[1.0, 0.0, 0.0], \
                     [0.0, np.cos(rollangle), -np.sin(rollangle)], \
                     [0.0, np.sin(rollangle), np.cos(rollangle)]])
        # Apply roll rotation around rotation centre
        xx=np.reshape(Xp0,(1,(m+1)*(n+1)),order='C')-xf
        yy=np.reshape(Yp0,(1,(m+1)*(n+1)),order='C')-yf
        zz=np.reshape(Zp0,(1,(m+1)*(n+1)),order='C')-zf
        dummy=Rx @ np.concatenate((xx,yy,zz),axis=0)
        # Assign rotated coordinates
        Xp0=np.reshape(dummy[0,:],(m+1,n+1),order='C')+xf
        Yp0=np.reshape(dummy[1,:],(m+1,n+1),order='C')+yf
        Zp0=np.reshape(dummy[2,:],(m+1,n+1),order='C')+zf
    # End if
    
    # Pitch the entire wing
    if rollpitchyaw[1] != 0.0 :
        # Rotation angle
        pitchangle=rollpitchyaw[1]
        # Rotation centre
        xf=rollpitchyaw_cent[0]
        yf=rollpitchyaw_cent[1]
        zf=rollpitchyaw_cent[2]
        # Rotation matrix
        Ry=np.array([[np.cos(pitchangle), 0.0, np.sin(pitchangle)], \
                     [0.0, 1.0, 0.0], \
                     [-np.sin(pitchangle), 0.0, np.cos(pitchangle)]])
        # Apply pitch rotation around rotation centre
        xx=np.reshape(Xp0,(1,(m+1)*(n+1)),order='C')-xf
        yy=np.reshape(Yp0,(1,(m+1)*(n+1)),order='C')-yf
        zz=np.reshape(Zp0,(1,(m+1)*(n+1)),order='C')-zf
        dummy=Ry @ np.concatenate((xx,yy,zz),axis=0)
        # Assign rotated coordinates
        Xp0=np.reshape(dummy[0,:],(m+1,n+1),order='C')+xf
        Yp0=np.reshape(dummy[1,:],(m+1,n+1),order='C')+yf
        Zp0=np.reshape(dummy[2,:],(m+1,n+1),order='C')+zf
    # End if

    # Yaw the entire wing
    if rollpitchyaw[2] != 0.0 :
        # Rotation angle
        yawangle=rollpitchyaw[2]
        # Rotation centre
        xf=rollpitchyaw_cent[0]
        yf=rollpitchyaw_cent[1]
        zf=rollpitchyaw_cent[2]
        # Rotation matrix
        Rz=np.array([[np.cos(yawangle), -np.sin(yawangle),0.0], \
                     [np.sin(yawangle), np.cos(yawangle),0.0], \
                     [0.0, 0.0, 1.0]])
        # Apply yaw rotation around rotation centre
        xx=np.reshape(Xp0,(1,(m+1)*(n+1)),order='C')-xf
        yy=np.reshape(Yp0,(1,(m+1)*(n+1)),order='C')-yf
        zz=np.reshape(Zp0,(1,(m+1)*(n+1)),order='C')-zf
        dummy=Rz @ np.concatenate((xx,yy,zz),axis=0)
        # Assign rotated coordinates
        Xp0=np.reshape(dummy[0,:],(m+1,n+1),order='C')+xf
        Yp0=np.reshape(dummy[1,:],(m+1,n+1),order='C')+yf
        Zp0=np.reshape(dummy[2,:],(m+1,n+1),order='C')+zf
    # End if
    
    # Non-dimensional spanwise control points coordinates
    yc=(yp[1:n+1]+yp[0:n])/2.0
    
    # Calculate control points, normal vectors and areas
    Xc0,Yc0,Zc0,nx0,ny0,nz0,tauxx0,tauxy0,tauxz0,tauyx0,tauyy0,tauyz0,s0,cpln0,sm0,tmx0,tmy0,tmz0,sn0,tnx0,tny0,tnz0=SDPM_control_normal_tangential(Xp0,Yp0,Zp0,dir_tau)
   
    # Reshape control point coordinate matrices into vectors
    Xc0all=np.reshape(Xc0,(m*n,1),order='C')
    Yc0all=np.reshape(Yc0,(m*n,1),order='C')
    Zc0all=np.reshape(Zc0,(m*n,1),order='C')
    # Reshape normal vector component matrices into vectors
    nx0all=np.reshape(nx0,(m*n,1),order='C')
    ny0all=np.reshape(ny0,(m*n,1),order='C')
    nz0all=np.reshape(nz0,(m*n,1),order='C')
    # Reshape panel area matrix into a vector
    s0all=np.reshape(s0,(m*n,1),order='C')

    # Create wake
    Xw0,Yw0,Zw0,mw,wake_dx=create_wake(nchords,Dxw,Xp0,Yp0,Zp0,cref)
    
    # Assign data to body struct array
    body['Xp0'][ibody]=Xp0
    body['Yp0'][ibody]=Yp0
    body['Zp0'][ibody]=Zp0
    body['Xc0'][ibody]=Xc0
    body['Yc0'][ibody]=Yc0
    body['Zc0'][ibody]=Zc0 
    body['Xc0all'][ibody]=Xc0all
    body['Yc0all'][ibody]=Yc0all
    body['Zc0all'][ibody]=Zc0all 
    body['Xw0'][ibody]=Xw0
    body['Yw0'][ibody]=Yw0
    body['Zw0'][ibody]=Zw0
    body['mw'][ibody]=mw
    body['wake_dx'][ibody]=wake_dx
    body['nx0'][ibody]=nx0
    body['ny0'][ibody]=ny0
    body['nz0'][ibody]=nz0
    body['tauxx0'][ibody]=tauxx0
    body['tauxy0'][ibody]=tauxy0
    body['tauxz0'][ibody]=tauxz0
    body['tauyx0'][ibody]=tauyx0
    body['tauyy0'][ibody]=tauyy0
    body['tauyz0'][ibody]=tauyz0
    body['nx0all'][ibody]=nx0all
    body['ny0all'][ibody]=ny0all
    body['nz0all'][ibody]=nz0all
    body['s0'][ibody]=s0
    body['s0all'][ibody]=s0all
    body['m'][ibody]=m
    body['n'][ibody]=n
    body['c0'][ibody]=trap['rootchord'][0]
    body['b'][ibody]=2*bhalf
    body['yc'][ibody]=yc
    body['AR'][ibody]=AR
    body['S'][ibody]=S
    body['name'][ibody]=name
    body['dir_tau'][ibody]=dir_tau

    return body

def create_wake(nchords,Dxw,Xp0,Yp0,Zp0,cref):
    # Create wake for the SDPM method, shed from the trailing edge of wing
    # whose coordinates are given in Xp0,Yp0,Zp0.
    # nchords: Length of wake in mean chords.
    # Dxw: Chordwise wake spacing parameter. Dxw=0: No wake, Dxw=-1: One 
    #        chordwise panel (steady case), Dxw>0: The chordwise wake spacing 
    #        will be Dxw  times the mean chordwise spacing of the corresponding 
    #        wing.
    # Xp0: (m+1)*(n+1) matrix of x-coordinates of wing panel vertices
    # Yp0: (m+1)*(n+1) matrix of y-coordinates of wing panel vertices
    # Zp0: (m+1)*(n+1) matrix of z-coordinates of wing panel vertices
    # cref: Reference chord length (usually mean aerodynamic chord)
    # Xw0: (mw+1)*(n+1) matrix of x-coordinates of wake panel vertices
    # Yw0: (mw+1)*(n+1) matrix of y-coordinates of wake panel vertices
    # Zw0: (mw+1)*(n+1) matrix of z-coordinates of wake panel vertices
    # mw: Number of chordwise wake panels
    # wake_dx: (mw+1)*1 vector of chordwise spacing of wake panels

    # Get dimensions of wing coordinate arrays
    m,n=Xp0.shape
    m=m-1       # Number of chordwise panels in body
    mhalf=m//2  # m is even for all wing-type bodies
    n=n-1       # Number of spanwise panels in body
    # Select case
    if Dxw == 0.0:    # This body sheds no wake
        mw=0          # Number of chordwise wake panels 
        Xw0=0.0       # x-coordinates of wake doublet panel vertices
        Yw0=0.0       # y-coordinates of wake doublet panel vertices
        Zw0=0.0       # z-coordinates of wake doublet panel vertices
        wake_dx=0.0   #  Wake chordwise spacing vector for all wake panels
    elif Dxw == -1.0: # One chordwise wake panel (steady case)
        mw=1                        # Number of chordwise wake panels 
        wake_dx=nchords*cref        # Chordwise length of wake
        Xw0=np.zeros((mw+1,n+1))       # Initialize x-coordinates of wake doublet panel vertices
        Xw0[0,:]=Xp0[m,:]           # Leading edge of wake doublet panels
        Xw0[1:mw+1,:]=Xw0[0,:]+wake_dx   # Trailing edge of wake doublet panels
        # Set y-coordinates of all spanwise wake doublet panels to corresponding
        # spanwise y-coordinates of wing panels
        Yw0=np.matlib.repmat(Yp0[m,:],mw+1,1)
        # Set z-coordinates of all spanwise wake doublet panels to corresponding
        # spanwise z-coordinates of wing panels
        Zw0=np.matlib.repmat(Zp0[m,:],mw+1,1)
    else:
        dx=Dxw*cref/mhalf                   # Chordwise length of wake doublet panels, use root chord to calculate them
        mw=int(np.floor(nchords*cref/dx))   # Number of chordwise wake panels
        wake_dx=np.ones((mw,1))*dx              # Wake spacing vector for all chordwise wake panels
        Xw0=np.zeros((mw+1,n+1))       # Initialize x-coordinates of wake doublet panel vertices
        Xw0[0,:]=Xp0[m,:]           # Leading edge of wake doublet panels
        Xw0[1:mw+1,:]=Xw0[0,:]+np.matlib.repmat(np.cumsum(wake_dx,axis=0),1,n+1)   # Trailing edge of wake doublet panels
        # Set y-coordinates of all spanwise wake doublet panels to corresponding
        # spanwise y-coordinates of wing panels
        Yw0=np.matlib.repmat(Yp0[m,:],mw+1,1)
        # Set z-coordinates of all spanwise wake doublet panels to corresponding
        # spanwise z-coordinates of wing panels
        Zw0=np.matlib.repmat(Zp0[m,:],mw+1,1)
    # End if
    return Xw0,Yw0,Zw0,mw,wake_dx
        
        
def SDPM_control_normal_tangential(Xp,Yp,Zp,dir_tau):   
    # Calculates the control points, normal unit vectors, tangent unit vectors,
    # surface areas, flatness, chordwise lengths, chordwise unit vectors,
    # spanwise lengths and spanwise unit vectors of the quadrilateral panels
    # whose endpoints are given in matrices Xp, Yp and Zp. The control points
    # lie on the mean of the x, y and z coordinates of each panel; this code is
    # therefore only compatiple with the Source and Doublet Panel Method. It is
    # not compatible with the Doublet Lattice Method or the Vortex Lattice
    # Method.
    #
    # Xp: x-coordinates of panel vertices arranged in a (m+1)*(n+1) matrix
    # Yp: y-coordinates of panel vertices arranged in a (m+1)*(n+1) matrix
    # Zp: z-coordinates of panel vertices arranged in a (m+1)*(n+1) matrix
    # dir: the direction in which the unit tangent vector (tauxx, tauxy, tauxz) has a zero component
    # Xc: x-coordinates of panel control points arranged in a m*n matrix
    # Yc: y-coordinates of panel control points arranged in a m*n matrix
    # Zc: z-coordinates of panel control points arranged in a m*n matrix
    # nx: x-coordinates of normal unit vectors arranged in a m*n matrix
    # ny: y-coordinates of normal unit vectors arranged in a m*n matrix
    # nz: z-coordinates of normal unit vectors arranged in a m*n matrix
    # tauxx: x-components of tangent unit vectors with zero y component
    #        arranged in a m*n matrix
    # tauxy: y-components of tangent unit vectors with zero y component
    #        arranged in a m*n matrix
    # tauxz: z-components of tangent unit vectors with zero y component
    #        arranged in a m*n matrix
    # tauyx: x-components of tangent unit vectors perpendicular to normal
    #        vector and first tangent vector, arranged in a m*n matrix
    # tauyy: y-components of tangent unit vectors perpendicular to normal
    #        vector and first tangent vector, arranged in a m*n matrix
    # tauyz: z-components of tangent unit vectors perpendicular to normal
    #        vector and first tangent vector, arranged in a m*n matrix
    # s: panel areas arranged in a m*n matrix
    # cpln: panel flatness parameter arranged in a m*n matrix. This is
    # calculated from the cross product of two normal unit vectors on each
    # panel evaluated from different vertices. If the two unit vectors are
    # colinear, their cross product is zero and the panel is flat. If the panel
    # is not flat, cpln is a measure of its flatness.
    # sm: m*n matrix of chordwise panel lengths
    # tmx: x-components of chordwise tangent unit vectors arranged in a m*n
    #      matrix
    # tmy: y-components of chordwise tangent unit vectors arranged in a m*n
    #      matrix
    # tmz: z-components of chordwise tangent unit vectors arranged in a m*n
    #      matrix
    # sn: m*n matrix of spanwise panel lengths
    # tnx: x-components of spanwise tangent unit vectors arranged in a m*n
    #      matrix
    # tny: y-components of spanwise tangent unit vectors arranged in a m*n
    #      matrix
    # tnz: z-components of spanwise tangent unit vectors arranged in a m*n
    #      matrix
    #
    # If the panels are not perfectly flat, the control points will not lie on
    # them. This is OK if the cpln is no too high. The code will probably work
    # for triangular panels but four vertices must still be defined, two of
    # them being identical.
    
    # Calculate m+1 and n+1
    mp1,np1=Xp.shape
    # Calculate m and n
    m=mp1-1
    n=np1-1
    # Set up all arrays. Xc, Yc and Zc are set up as matrices because I 
    # can't concatenate arrays properly yet. They are converted to arrays before
    # being returned
    Xc=np.matrix(np.zeros((m, n)))
    Yc=np.matrix(np.zeros((m, n)))
    Zc=np.matrix(np.zeros((m, n)))
    nx=np.zeros((m, n))
    ny=np.zeros((m, n))
    nz=np.zeros((m, n))
    tauxx=np.zeros((m, n))
    tauxy=np.zeros((m, n))
    tauxz=np.zeros((m, n))
    tauyx=np.zeros((m, n))
    tauyy=np.zeros((m, n))
    tauyz=np.zeros((m, n))
    s=np.zeros((m, n))
    cpln=np.zeros((m, n))
    # Calculate control points, surfaces, normal and tangent vectors
    for i in range(0,m): 
        for j in range(0,n): 
            # Calculate control points
            x=np.array([Xp[i,j], Xp[i,j+1], Xp[i+1,j+1], Xp[i+1,j]])
            y=np.array([Yp[i,j], Yp[i,j+1], Yp[i+1,j+1], Yp[i+1,j]])
            z=np.array([Zp[i,j], Zp[i,j+1], Zp[i+1,j+1], Zp[i+1,j]])
            Xc[i,j]=np.mean(x)
            Yc[i,j]=np.mean(y)
            Zc[i,j]=np.mean(z)
            # Calculate the normal vectors and panel areas
            N,ss,cplnn=normalarea(x,y,z)
            nx[i,j]  = N[0]
            ny[i,j]  = N[1]
            nz[i,j]  = N[2]
            s[i,j] = ss
            cpln[i,j] = cplnn
            
            if dir_tau == 1:
                # Calculate tangential vector with zero x-component
                tauxx[i,j]=0.0                              # x-component of vector
                tauxy[i,j]=1/np.sqrt(1+(N[1]/N[2])**2)      # y-component of vector
                tauxz[i,j]=-N[1]/N[2]*tauxy[i,j]            # z-component of vector
            elif dir_tau == 2:
                # Calculate tangential vector with zero y-component
                tauxx[i,j]=1/np.sqrt(1+(N[0]/N[2])**2)      # x-component of vector
                tauxy[i,j]=0.0                              # y-component of vector
                tauxz[i,j]=-N[0]/N[2]*tauxx[i,j]            # z-component of vector
            elif dir_tau == 3:
                # Calculate tangential vector with zero z-component
                tauxx[i,j]=1/np.sqrt(1+(N[0]/N[1])**2)      # x-component of vector
                tauxy[i,j]=-N[0]/N[1]*tauxx[i,j]            # y-component of vector
                tauxz[i,j]=0.0                              # z-component of vector
            # End if     

            # Calculate a second tangential vector, perpendicular to the normal
            # and first tangential vector from equation 5.185
            tauy_xyz=np.cross(N,np.array([tauxx[i,j], tauxy[i,j], tauxz[i,j]]))
            tauyx[i,j]=tauy_xyz[0]     # x-component of vector
            tauyy[i,j]=tauy_xyz[1]     # y-component of vector
            tauyz[i,j]=tauy_xyz[2]     # z-component of vector
        # End if
    # End if
    
    # Now calculate tangent vectors and lengths necessary for calculation of surface velocities
    if m > 1:        
        # Calculate chordwise lengths of panels. They are also converted to arrays from matrices
        dXm=np.array(np.concatenate((Xc[1,:]-Xc[0,:],Xc[2:m,:]-Xc[0:m-2,:],Xc[m-1,:]-Xc[m-2,:]),axis=0))  # In x direction
        dYm=np.array(np.concatenate((Yc[1,:]-Yc[0,:],Yc[2:m,:]-Yc[0:m-2,:],Yc[m-1,:]-Yc[m-2,:]),axis=0))  # In y direction
        dZm=np.array(np.concatenate((Zc[1,:]-Zc[0,:],Zc[2:m,:]-Zc[0:m-2,:],Zc[m-1,:]-Zc[m-2,:]),axis=0))  # In z direction
        sm=np.sqrt(np.square(dXm)+np.square(dYm)+np.square(dZm))      # Total length
        # Calculate chordwise tangent unit vectors from equations 5.191
        tmx=np.divide(dXm,sm)    # x-component
        tmy=np.divide(dYm,sm)    # y-component
        tmz=np.divide(dZm,sm)    # z-component
    else:
        # It is impossible to calculate derivatives in the chordwise
        # direction
        sm=np.ones((m,n))
        tmx=np.zeros((m,n))
        tmy=np.zeros((m,n))
        tmz=np.zeros((m,n))        
    # End if
        
    # Calculate spanwise lengths of panels.  They are also converted to arrays from matrices
    if n > 2:
        dXn=np.array(np.concatenate((Xc[:,1]-Xc[:,0],Xc[:,2:n]-Xc[:,0:n-2],Xc[:,n-1]-Xc[:,n-2]),axis=1))  # In x direction
        dYn=np.array(np.concatenate((Yc[:,1]-Yc[:,0],Yc[:,2:n]-Yc[:,0:n-2],Yc[:,n-1]-Yc[:,n-2]),axis=1))  # In x direction
        dZn=np.array(np.concatenate((Zc[:,1]-Zc[:,0],Zc[:,2:n]-Zc[:,0:n-2],Zc[:,n-1]-Zc[:,n-2]),axis=1))  # In x direction
        sn=np.sqrt(np.square(dXn)+np.square(dYn)+np.square(dZn))      # Total length
        # Calculate spanwise tangent unit vectors from equations 5.193
        tnx=np.divide(dXn,sn)    # x-component
        tny=np.divide(dYn,sn)    # x-component
        tnz=np.divide(dZn,sn)    # x-component
    elif n == 2:
        dXn=np.array(np.concatenate((Xc[:,1]-Xc[:,0],Xc[:,n-1]-Xc[:,n-2]),axis=1))  # In x direction
        dYn=np.array(np.concatenate((Yc[:,1]-Yc[:,0],Yc[:,n-1]-Yc[:,n-2]),axis=1))  # In x direction
        dZn=np.array(np.concatenate((Zc[:,1]-Zc[:,0],Zc[:,n-1]-Zc[:,n-2]),axis=1))  # In x direction
        sn=np.sqrt(np.square(dXn)+np.square(dYn)+np.square(dZn))      # Total length
        # Calculate spanwise tangent unit vectors from equations 5.193
        tnx=np.divide(dXn,sn)    # x-component
        tny=np.divide(dYn,sn)    # x-component
        tnz=np.divide(dZn,sn)    # x-component
    else:
        # It is impossible to calculate derivatives in the spanwise
        # direction
        sn=np.ones((m,n))
        tnx=np.zeros((m,n))
        tny=np.zeros((m,n))
        tnz=np.zeros((m,n))
    # End if    
    
    Xc=np.array(Xc)
    Yc=np.array(Yc)
    Zc=np.array(Zc)

    return Xc,Yc,Zc,nx,ny,nz,tauxx,tauxy,tauxz,tauyx,tauyy,tauyz,s,cpln,sm,tmx,tmy,tmz,sn,tnx,tny,tnz

def makewingtips(body,ibody,ibody0,mirroredwing):
    # Creates SDPM grids for wingtips for the wing stored in body(ibody).
    # body: struct array describing the bodies present in the flow
    # ibody: index of body for which to create wingtips
    # ibody0: index of body after which to store wingtips
    # mirroredwing: determines if the body(ibody) is a left-sided, right-sided
    #       or two-sided wing
    # IMPORTANT: makewingtips stors the winglets in positions ibody0+1 (and
    # ibody0+2 if mirroredwing=2) of the body struct array. If any other data 
    # is stored in these positions it will be over-written. 
    
    # Get number of chordwise panels
    m=body['m'][ibody]//2
    # Set number of heightwise panels
    n=2
    # Determine direction of taux tangential vector
    if body['dir_tau'][ibody] == 2:
        dir_tau=3
    elif body['dir_tau'][ibody] == 3:
        dir_tau=1
    # End if
    
    if mirroredwing == 2:
        # Create left wingtip
        # Acquire leftmost panel vertex coordinates 
        X=body['Xp0'][ibody][:,0]
        Y=body['Yp0'][ibody][:,0]
        Z=body['Zp0'][ibody][:,0]
        # Calculate tip chord
        c0=np.sqrt((X[0]-X[m])**2+(Y[0]-Y[m])**2+(Z[0]-Z[m])**2)
        # Create SDPM grid for left wingtip
        Xp0=np.zeros((m+1,3))
        Yp0=np.zeros((m+1,3))
        Zp0=np.zeros((m+1,3))
        Xp0[:,0]=np.flipud(X[0:m+1])
        Xp0[:,1]=(np.flipud(X[0:m+1])+X[m:2*m+1])/2.0
        Xp0[:,2]=X[m:2*m+1]
        Yp0[:,0]=np.flipud(Y[0:m+1])
        Yp0[:,1]=(np.flipud(Y[0:m+1])+Y[m:2*m+1])/2.0
        Yp0[:,2]=Y[m:2*m+1]
        Zp0[:,0]=np.flipud(Z[0:m+1])
        Zp0[:,1]=(np.flipud(Z[0:m+1])+Z[m:2*m+1])/2.0
        Zp0[:,2]=Z[m:2*m+1]
        # Calculate control points, normal vectors and areas
        Xc0,Yc0,Zc0,nx0,ny0,nz0,tauxx0,tauxy0,tauxz0,tauyx0,tauyy0,tauyz0,s0,cpln0,sm0,tmx0,tmy0,tmz0,sn0,tnx0,tny0,tnz0=SDPM_control_normal_tangential(Xp0,Yp0,Zp0,dir_tau)
        
        # Reshape control point coordinate matrices into vectors
        Xc0all=np.reshape(Xc0,(m*n,1),order='C')
        Yc0all=np.reshape(Yc0,(m*n,1),order='C')
        Zc0all=np.reshape(Zc0,(m*n,1),order='C')
        # Reshape normal vector component matrices into vectors
        nx0all=np.reshape(nx0,(m*n,1),order='C')
        ny0all=np.reshape(ny0,(m*n,1),order='C')
        nz0all=np.reshape(nz0,(m*n,1),order='C')
        # Reshape panel area matrix into a vector
        s0all=np.reshape(s0,(m*n,1),order='C')
        
        # Assign data to body struct array
        body['Xp0'][ibody0+1]=Xp0
        body['Yp0'][ibody0+1]=Yp0
        body['Zp0'][ibody0+1]=Zp0
        body['Xc0'][ibody0+1]=Xc0
        body['Yc0'][ibody0+1]=Yc0
        body['Zc0'][ibody0+1]=Zc0 
        body['Xc0all'][ibody0+1]=Xc0all
        body['Yc0all'][ibody0+1]=Yc0all
        body['Zc0all'][ibody0+1]=Zc0all 
        body['mw'][ibody0+1]=0
        body['nx0'][ibody0+1]=nx0
        body['ny0'][ibody0+1]=ny0
        body['nz0'][ibody0+1]=nz0
        body['tauxx0'][ibody0+1]=tauxx0
        body['tauxy0'][ibody0+1]=tauxy0
        body['tauxz0'][ibody0+1]=tauxz0
        body['tauyx0'][ibody0+1]=tauyx0
        body['tauyy0'][ibody0+1]=tauyy0
        body['tauyz0'][ibody0+1]=tauyz0
        body['nx0all'][ibody0+1]=nx0all
        body['ny0all'][ibody0+1]=ny0all
        body['nz0all'][ibody0+1]=nz0all
        body['s0'][ibody0+1]=s0
        body['s0all'][ibody0+1]=s0all
        body['m'][ibody0+1]=m
        body['n'][ibody0+1]=n
        body['c0'][ibody0+1]=c0
        body['name'][ibody0+1]='left_wingtip'
        body['dir_tau'][ibody0+1]=dir_tau

        # Create right wingtip
        # Acquire rightmost panel vertex coordinates 
        X=body['Xp0'][ibody][:,body['n'][ibody]]
        Y=body['Yp0'][ibody][:,body['n'][ibody]]
        Z=body['Zp0'][ibody][:,body['n'][ibody]]
        # Calculate tip chord
        c0=np.sqrt((X[0]-X[m])**2+(Y[0]-Y[m])**2+(Z[0]-Z[m])**2)
        # Create SDPM grid for left wingtip
        Xp0=np.zeros((m+1,3))
        Yp0=np.zeros((m+1,3))
        Zp0=np.zeros((m+1,3))
        Xp0[:,0]=X[m:2*m+1]
        Xp0[:,1]=(np.flipud(X[0:m+1])+X[m:2*m+1])/2.0
        Xp0[:,2]=np.flipud(X[0:m+1])
        Yp0[:,0]=Y[m:2*m+1]
        Yp0[:,1]=(np.flipud(Y[0:m+1])+Y[m:2*m+1])/2.0
        Yp0[:,2]=np.flipud(Y[0:m+1])
        Zp0[:,0]=Z[m:2*m+1]
        Zp0[:,1]=(np.flipud(Z[0:m+1])+Z[m:2*m+1])/2.0
        Zp0[:,2]=np.flipud(Z[0:m+1])
        # Calculate control points, normal vectors and areas
        Xc0,Yc0,Zc0,nx0,ny0,nz0,tauxx0,tauxy0,tauxz0,tauyx0,tauyy0,tauyz0,s0,cpln0,sm0,tmx0,tmy0,tmz0,sn0,tnx0,tny0,tnz0=SDPM_control_normal_tangential(Xp0,Yp0,Zp0,dir_tau)
        
        # Reshape control point coordinate matrices into vectors
        Xc0all=np.reshape(Xc0,(m*n,1),order='C')
        Yc0all=np.reshape(Yc0,(m*n,1),order='C')
        Zc0all=np.reshape(Zc0,(m*n,1),order='C')
        # Reshape normal vector component matrices into vectors
        nx0all=np.reshape(nx0,(m*n,1),order='C')
        ny0all=np.reshape(ny0,(m*n,1),order='C')
        nz0all=np.reshape(nz0,(m*n,1),order='C')
        # Reshape panel area matrix into a vector
        s0all=np.reshape(s0,(m*n,1),order='C')
        
        # Assign data to body struct array
        body['Xp0'][ibody0+2]=Xp0
        body['Yp0'][ibody0+2]=Yp0
        body['Zp0'][ibody0+2]=Zp0
        body['Xc0'][ibody0+2]=Xc0
        body['Yc0'][ibody0+2]=Yc0
        body['Zc0'][ibody0+2]=Zc0 
        body['Xc0all'][ibody0+2]=Xc0all
        body['Yc0all'][ibody0+2]=Yc0all
        body['Zc0all'][ibody0+2]=Zc0all 
        body['mw'][ibody0+2]=0
        body['nx0'][ibody0+2]=nx0
        body['ny0'][ibody0+2]=ny0
        body['nz0'][ibody0+2]=nz0
        body['tauxx0'][ibody0+2]=tauxx0
        body['tauxy0'][ibody0+2]=tauxy0
        body['tauxz0'][ibody0+2]=tauxz0
        body['tauyx0'][ibody0+2]=tauyx0
        body['tauyy0'][ibody0+2]=tauyy0
        body['tauyz0'][ibody0+2]=tauyz0
        body['nx0all'][ibody0+2]=nx0all
        body['ny0all'][ibody0+2]=ny0all
        body['nz0all'][ibody0+2]=nz0all
        body['s0'][ibody0+2]=s0
        body['s0all'][ibody0+2]=s0all
        body['m'][ibody0+2]=m
        body['n'][ibody0+2]=n
        body['c0'][ibody0+2]=c0
        body['name'][ibody0+2]='right_wingtip'
        body['dir_tau'][ibody0+2]=dir_tau     
    elif mirroredwing == -1:
        # Create left wingtip
        # Acquire leftmost panel vertex coordinates 
        X=body['Xp0'][ibody][:,0]
        Y=body['Yp0'][ibody][:,0]
        Z=body['Zp0'][ibody][:,0]
        # Calculate tip chord
        c0=np.sqrt((X[0]-X[m])**2+(Y[0]-Y[m])**2+(Z[0]-Z[m])**2)
        # Create SDPM grid for left wingtip
        Xp0=np.zeros((m+1,3))
        Yp0=np.zeros((m+1,3))
        Zp0=np.zeros((m+1,3))
        Xp0[:,0]=np.flipud(X[0:m+1])
        Xp0[:,1]=(np.flipud(X[0:m+1])+X[m:2*m+1])/2.0
        Xp0[:,2]=X[m:2*m+1]
        Yp0[:,0]=np.flipud(Y[0:m+1])
        Yp0[:,1]=(np.flipud(Y[0:m+1])+Y[m:2*m+1])/2.0
        Yp0[:,2]=Y[m:2*m+1]
        Zp0[:,0]=np.flipud(Z[0:m+1])
        Zp0[:,1]=(np.flipud(Z[0:m+1])+Z[m:2*m+1])/2.0
        Zp0[:,2]=Z[m:2*m+1]
        # Calculate control points, normal vectors and areas
        Xc0,Yc0,Zc0,nx0,ny0,nz0,tauxx0,tauxy0,tauxz0,tauyx0,tauyy0,tauyz0,s0,cpln0,sm0,tmx0,tmy0,tmz0,sn0,tnx0,tny0,tnz0=SDPM_control_normal_tangential(Xp0,Yp0,Zp0,dir_tau)
        
        # Reshape control point coordinate matrices into vectors
        Xc0all=np.reshape(Xc0,(m*n,1),order='C')
        Yc0all=np.reshape(Yc0,(m*n,1),order='C')
        Zc0all=np.reshape(Zc0,(m*n,1),order='C')
        # Reshape normal vector component matrices into vectors
        nx0all=np.reshape(nx0,(m*n,1),order='C')
        ny0all=np.reshape(ny0,(m*n,1),order='C')
        nz0all=np.reshape(nz0,(m*n,1),order='C')
        # Reshape panel area matrix into a vector
        s0all=np.reshape(s0,(m*n,1),order='C')
        
        # Assign data to body struct array
        body['Xp0'][ibody0+1]=Xp0
        body['Yp0'][ibody0+1]=Yp0
        body['Zp0'][ibody0+1]=Zp0
        body['Xc0'][ibody0+1]=Xc0
        body['Yc0'][ibody0+1]=Yc0
        body['Zc0'][ibody0+1]=Zc0 
        body['Xc0all'][ibody0+1]=Xc0all
        body['Yc0all'][ibody0+1]=Yc0all
        body['Zc0all'][ibody0+1]=Zc0all 
        body['mw'][ibody0+1]=0
        body['nx0'][ibody0+1]=nx0
        body['ny0'][ibody0+1]=ny0
        body['nz0'][ibody0+1]=nz0
        body['tauxx0'][ibody0+1]=tauxx0
        body['tauxy0'][ibody0+1]=tauxy0
        body['tauxz0'][ibody0+1]=tauxz0
        body['tauyx0'][ibody0+1]=tauyx0
        body['tauyy0'][ibody0+1]=tauyy0
        body['tauyz0'][ibody0+1]=tauyz0
        body['nx0all'][ibody0+1]=nx0all
        body['ny0all'][ibody0+1]=ny0all
        body['nz0all'][ibody0+1]=nz0all
        body['s0'][ibody0+1]=s0
        body['s0all'][ibody0+1]=s0all
        body['m'][ibody0+1]=m
        body['n'][ibody0+1]=n
        body['c0'][ibody0+1]=c0
        body['name'][ibody0+1]='left_wingtip'
        body['dir_tau'][ibody0+1]=dir_tau
    elif mirroredwing == 1:
        # Create right wingtip
        # Acquire rightmost panel vertex coordinates 
        X=body['Xp0'][ibody][:,body['n'][ibody]]
        Y=body['Yp0'][ibody][:,body['n'][ibody]]
        Z=body['Zp0'][ibody][:,body['n'][ibody]]
        # Calculate tip chord
        c0=np.sqrt((X[0]-X[m])**2+(Y[0]-Y[m])**2+(Z[0]-Z[m])**2)
        # Create SDPM grid for left wingtip
        Xp0=np.zeros((m+1,3))
        Yp0=np.zeros((m+1,3))
        Zp0=np.zeros((m+1,3))
        Xp0[:,0]=X[m:2*m+1]
        Xp0[:,1]=(np.flipud(X[0:m+1])+X[m:2*m+1])/2.0
        Xp0[:,2]=np.flipud(X[0:m+1])
        Yp0[:,0]=Y[m:2*m+1]
        Yp0[:,1]=(np.flipud(Y[0:m+1])+Y[m:2*m+1])/2.0
        Yp0[:,2]=np.flipud(Y[0:m+1])
        Zp0[:,0]=Z[m:2*m+1]
        Zp0[:,1]=(np.flipud(Z[0:m+1])+Z[m:2*m+1])/2.0
        Zp0[:,2]=np.flipud(Z[0:m+1])
        # Calculate control points, normal vectors and areas
        Xc0,Yc0,Zc0,nx0,ny0,nz0,tauxx0,tauxy0,tauxz0,tauyx0,tauyy0,tauyz0,s0,cpln0,sm0,tmx0,tmy0,tmz0,sn0,tnx0,tny0,tnz0=SDPM_control_normal_tangential(Xp0,Yp0,Zp0,dir_tau)
        
        # Reshape control point coordinate matrices into vectors
        Xc0all=np.reshape(Xc0,(m*n,1),order='C')
        Yc0all=np.reshape(Yc0,(m*n,1),order='C')
        Zc0all=np.reshape(Zc0,(m*n,1),order='C')
        # Reshape normal vector component matrices into vectors
        nx0all=np.reshape(nx0,(m*n,1),order='C')
        ny0all=np.reshape(ny0,(m*n,1),order='C')
        nz0all=np.reshape(nz0,(m*n,1),order='C')
        # Reshape panel area matrix into a vector
        s0all=np.reshape(s0,(m*n,1),order='C')
        
        # Assign data to body struct array
        body['Xp0'][ibody0+1]=Xp0
        body['Yp0'][ibody0+1]=Yp0
        body['Zp0'][ibody0+1]=Zp0
        body['Xc0'][ibody0+1]=Xc0
        body['Yc0'][ibody0+1]=Yc0
        body['Zc0'][ibody0+1]=Zc0 
        body['Xc0all'][ibody0+1]=Xc0all
        body['Yc0all'][ibody0+1]=Yc0all
        body['Zc0all'][ibody0+1]=Zc0all 
        body['mw'][ibody0+1]=0
        body['nx0'][ibody0+1]=nx0
        body['ny0'][ibody0+1]=ny0
        body['nz0'][ibody0+1]=nz0
        body['tauxx0'][ibody0+1]=tauxx0
        body['tauxy0'][ibody0+1]=tauxy0
        body['tauxz0'][ibody0+1]=tauxz0
        body['tauyx0'][ibody0+1]=tauyx0
        body['tauyy0'][ibody0+1]=tauyy0
        body['tauyz0'][ibody0+1]=tauyz0
        body['nx0all'][ibody0+1]=nx0all
        body['ny0all'][ibody0+1]=ny0all
        body['nz0all'][ibody0+1]=nz0all
        body['s0'][ibody0+1]=s0
        body['s0all'][ibody0+1]=s0all
        body['m'][ibody0+1]=m
        body['n'][ibody0+1]=n
        body['c0'][ibody0+1]=c0
        body['name'][ibody0+1]='right_wingtip'
        body['dir_tau'][ibody0+1]=dir_tau     
    # End if
    return body

def SDPMgeometry_fuselage(body,ibody,m,nlengths,Dxw,n,nmin,linlength,xfus,fuselage_sections,nosexyz,name):
    # Calculates the coordinates of fuselage panel vertices for the Source and
    # Doublet Panel method from trapezoidal section data.
    # This function is part of the SDPMflut Matlab distribution.
    # m:  Number of lengthwise panels (nose to tail).
    # nlengths: Length of wake in fuselage lengths
    # Dxw: Chordwise wake spacing parameter. Dxw=0: No wake, Dxw=-1: One 
    #        chordwise panel (steady case), Dxw>0: The chordwise wake spacing 
    #        will be Dxw  times the mean chordwise spacing of the corresponding 
    #        wing.
    # n: Number of heightwise panels (bottom to top)
    # nmin: Minimum number of panels per trapezoidal section
    # linlength: if linlength=1 the lengthsise panel distribution is linear,
    #           otherwise it is nonlinear
    # xfus: x-coordinates of trapezoidal sections
    # fuselage_sections: struct array containing the description of the
    #           sections of the fuselage.
    # Xp0: (m+1)*(n+1) matrix of x-coordinates of panel vertices
    # Yp0: (m+1)*(n+1) matrix of y-coordinates of panel vertices
    # Zp0: (m+1)*(n+1) matrix of z-coordinates of panel vertices
    # nosexyz: 1*3 vector of x,y,z coordinates of distance of nose leading edge 
    #       from (0,0,0)
    # name: string of name of body
    
    # Number of trapezoidal sections
    num_sec=len(xfus)-1
    # Total length of fuselage
    lfus=xfus[num_sec]-xfus[0]
    # Calculate lengths of trapezoidal sections
    trap_length=xfus[1:num_sec+1]-xfus[0:num_sec]
    # Cumulative length of trapezoidal sections
    lsum=np.block([0.0, np.cumsum(trap_length)])

    # Distribute lengthwise panels over the trapezoidal sections
    n_sec=np.round(trap_length/lfus*n) # Calculate number of lengthwise points per section
    n_sec=n_sec.astype(int)
    # Ensure that none of the sections have fewer than nmin spanwise panels
    iko=np.argwhere(n_sec <= nmin)
    if len(iko) > 0:
        n_sec[iko]=nmin
    # End if
    # Update lengthwise panel number
    n=np.sum(n_sec)

    # Initialize grid
    Xp0=np.zeros((m+1,n+1))
    Yp0=np.zeros((m+1,n+1))
    Zp0=np.zeros((m+1,n+1))
    
    # Heightwise grid in angle
    theta=np.linspace(-np.pi/2,np.pi/2,num=m+1) 
    # Resample all heightwise data at these angles
    Yres=np.zeros((m+1,num_sec+1))
    Zres=np.zeros((m+1,num_sec+1))
    for i in range (0,num_sec+1):
        # Get y and z coordinates of current fuselage section
        Y=fuselage_sections['y'][i][0]
        Z=fuselage_sections['z'][i][0]
        if len(Y)==0 or all(Y == 0):
            Yres[:,i]=Y[0]*np.ones(m+1)
            Zres[:,i]=Z[0]*np.ones(m+1)
        else:
            # Calculate midpoint of z coordinates
            Zmid=(Z[-1]-Z[0])/2.0+Z[0]
            # Transform to polar coordinates
            th=np.arctan2(Z-Zmid,Y)
            r=np.sqrt(Y**2+(Z-Zmid)**2)
            # Interpolate onto angles theta
            cs = CubicSpline(th, r)
            # Calculate airfoil shape at xp
            rint=cs(theta)
            # Transform to Cartesian coordinates
            Yint=rint*np.cos(theta)
            Zint=rint*np.sin(theta)
            # Yint must be exactly equal to zero at top and bottom
            Yint[0]=0.0
            Yint[m]=0.0
            # Store resampled y and z coordinates
            Yres[:,i]=Yint;
            Zres[:,i]=Zint+Zmid;
        # End if
    # End for
    
    # Cumulative sum of panels
    nsum=np.block([0, np.cumsum(n_sec)])
    # Cycle over the trapezoidal sections
    for i in range (0,num_sec):
        # Set up lengthwise grid
        if linlength == 1:
            xp=np.linspace(0,1,n_sec[i]+1)
        else:
            xp=(1.0+np.cos(np.linspace(np.pi,0.0,n_sec[i]+1)))/2.0
        # End if
        # Calculate y-coordinates of panel vertices (add xfus(1) to start at the right place)
        Xp0[:,nsum[i]:nsum[i+1]+1]=np.matlib.repmat(xp,m+1,1)*trap_length[i]+lsum[i]+xfus[0]   
        # Interpolate linearly the y coordinates
        slope=(Yres[:,i+1]-Yres[:,i])/(xfus[i+1]-xfus[i])
        Yp0[:,nsum[i]:nsum[i+1]+1]=np.matlib.repmat(np.reshape(slope,(m+1,1)),1,nsum[i+1]+1-nsum[i])* \
            np.matlib.repmat(Xp0[1,nsum[i]:nsum[i+1]+1]-xfus[i],m+1,1)+np.matlib.repmat(np.reshape(Yres[:,i],(m+1,1)),1,nsum[i+1]+1-nsum[i])
        # Interpolate linearly the x coordinates
        slope=(Zres[:,i+1]-Zres[:,i])/(xfus[i+1]-xfus[i])
        Zp0[:,nsum[i]:nsum[i+1]+1]=np.matlib.repmat(np.reshape(slope,(m+1,1)),1,nsum[i+1]+1-nsum[i])* \
            np.matlib.repmat(Xp0[1,nsum[i]:nsum[i+1]+1]-xfus[i],m+1,1)+np.matlib.repmat(np.reshape(Zres[:,i],(m+1,1)),1,nsum[i+1]+1-nsum[i])
    # End for
    
    # Move entire fuselage to required leading edge position
    Xp0=Xp0+nosexyz[0]
    Yp0=Yp0+nosexyz[1]
    Zp0=Zp0+nosexyz[2]
    
    # For fuselages dir_tau=1
    dir_tau=1
    # Calculate control points, normal vectors and areas
    Xc0,Yc0,Zc0,nx0,ny0,nz0,tauxx0,tauxy0,tauxz0,tauyx0,tauyy0,tauyz0,s0,cpln0,sm0,tmx0,tmy0,tmz0,sn0,tnx0,tny0,tnz0=SDPM_control_normal_tangential(Xp0,Yp0,Zp0,dir_tau)

    # Reshape control point coordinate matrices into vectors
    Xc0all=np.reshape(Xc0,(m*n,1),order='C')
    Yc0all=np.reshape(Yc0,(m*n,1),order='C')
    Zc0all=np.reshape(Zc0,(m*n,1),order='C')
    # Reshape normal vector component matrices into vectors
    nx0all=np.reshape(nx0,(m*n,1),order='C')
    ny0all=np.reshape(ny0,(m*n,1),order='C')
    nz0all=np.reshape(nz0,(m*n,1),order='C')
    # Reshape panel area matrix into a vector
    s0all=np.reshape(s0,(m*n,1),order='C')

    # Select case
    if Dxw == 0.0:    # This body sheds no wake
        mw=0          # Number of chordwise wake panels 
        Xw0=0.0       # x-coordinates of wake doublet panel vertices
        Yw0=0.0       # y-coordinates of wake doublet panel vertices
        Zw0=0.0       # z-coordinates of wake doublet panel vertices
        wake_dx=0.0   #  Wake chordwise spacing vector for all wake panels
    elif Dxw == -1.0: # One chordwise wake panel (steady case)
        mw=1                        # Number of chordwise wake panels 
        wake_dx=nlengths*lfus        # Chordwise length of wake
        Xw0=np.zeros((m+1,2))       # Initialize x-coordinates of wake doublet panel vertices
        Xw0[:,0]=Xp0[:,n]           # Leading edge of wake doublet panels
        Xw0[:,1]=Xw0[:,0]+wake_dx   # Trailing edge of wake doublet panels
        # Set y-coordinates of all spanwise wake doublet panels to corresponding
        # spanwise y-coordinates of wing panels
        Yw0=np.matlib.repmat(np.reshape(Yp0[:,n],(m+1,1)),1,mw+1)
        # Set z-coordinates of all spanwise wake doublet panels to corresponding
        # spanwise z-coordinates of wing panels
        Zw0=np.matlib.repmat(np.reshape(Zp0[:,n],(m+1,1)),1,mw+1)
    else:
        dx=Dxw*lfus/n                   # Chordwise length of wake doublet panels, use root chord to calculate them
        mw=int(np.floor(nlengths*lfus/dx))   # Number of chordwise wake panels
        wake_dx=np.ones((mw,1))*dx              # Wake spacing vector for all chordwise wake panels
        Xw0=np.zeros((m+1,mw+1))       # Initialize x-coordinates of wake doublet panel vertices
        Xw0[:,0]=Xp0[:,n]              # Leading edge of wake doublet panels
        Xw0[:,1:mw+1]=np.matlib.repmat(np.reshape(Xw0[:,0],(m+1,1)),1,mw)+np.matlib.repmat(np.cumsum(wake_dx),m+1,1)   # Trailing edge of wake doublet panels
        # Set y-coordinates of all spanwise wake doublet panels to corresponding
        # spanwise y-coordinates of wing panels
        Yw0=np.matlib.repmat(np.reshape(Yp0[:,n],(m+1,1)),1,mw+1)
        # Set z-coordinates of all spanwise wake doublet panels to corresponding
        # spanwise z-coordinates of wing panels
        Zw0=np.matlib.repmat(np.reshape(Zp0[:,n],(m+1,1)),1,mw+1)
    # End if

    # Assign data to body struct array
    body['Xp0'][ibody]=Xp0
    body['Yp0'][ibody]=Yp0
    body['Zp0'][ibody]=Zp0
    body['Xc0'][ibody]=Xc0
    body['Yc0'][ibody]=Yc0
    body['Zc0'][ibody]=Zc0 
    body['Xc0all'][ibody]=Xc0all
    body['Yc0all'][ibody]=Yc0all
    body['Zc0all'][ibody]=Zc0all 
    body['Xw0'][ibody]=Xw0
    body['Yw0'][ibody]=Yw0
    body['Zw0'][ibody]=Zw0
    body['mw'][ibody]=mw
    body['wake_dx'][ibody]=wake_dx
    body['nx0'][ibody]=nx0
    body['ny0'][ibody]=ny0
    body['nz0'][ibody]=nz0
    body['tauxx0'][ibody]=tauxx0
    body['tauxy0'][ibody]=tauxy0
    body['tauxz0'][ibody]=tauxz0
    body['tauyx0'][ibody]=tauyx0
    body['tauyy0'][ibody]=tauyy0
    body['tauyz0'][ibody]=tauyz0
    body['nx0all'][ibody]=nx0all
    body['ny0all'][ibody]=ny0all
    body['nz0all'][ibody]=nz0all
    body['s0'][ibody]=s0
    body['s0all'][ibody]=s0all
    body['m'][ibody]=m
    body['n'][ibody]=n
    body['c0'][ibody]=0.0
    body['b'][ibody]=lfus
    body['yc'][ibody]=0.0
    body['AR'][ibody]=0.0
    body['S'][ibody]=0.0
    body['name'][ibody]=name
    body['dir_tau'][ibody]=dir_tau

    return body

def in_triangle(xl,yl,xat0,yat0,tol):
    # Determines if point xat0,yat0 lies inside triangle defined by vertices in
    # xl,yl, allowing a tolerance tol.
    # Based on code in https://www.mathworks.com/matlabcentral/answers/277984-check-points-inside-triangle-or-on-edge-with-example
    # xl: 1x3 array of x coordinates of triangle vertices
    # yl: 1x3 array of y coordinates of triangle vertices
    # zl: 1x3 array of z coordinates of triangle vertices
    # xat0, yat0: x and y coordinates of point
    # tol: small positive real number. If point lies to within -tol of
    #       triangle, it is assumed to lie inside.
    # in: Logical 1 or 0. 1 if inside, 0 if outside.
    # res: mminimum of the three conditions
    
    # Assemble the three vertices of the triangle in vectors r1, r2, r3
    r1=np.array([xl[0], yl[0]])
    r2=np.array([xl[1], yl[1]])
    r3=np.array([xl[2], yl[2]])
    # Assemble point xat0, yat0 in vector r0
    r0=np.array([xat0, yat0])
    # Calculate distance between vertices
    r12=r1-r2
    r23=r2-r3
    r31=r3-r1
    # Initialize condition
    cond=np.zeros(3)
    # Condition for lying on or inside segment 23
    cond[0]=np.linalg.det(np.block([[r31], [r23]]))*np.linalg.det(np.block([[r3-r0], [r23]]))
    # Condition for lying on or inside segment 31
    cond[1]=np.linalg.det(np.block([[r12], [r31]]))*np.linalg.det(np.block([[r1-r0], [r31]]))
    # Condition for lying on or inside segment 12
    cond[2]=np.linalg.det(np.block([[r23], [r12]]))*np.linalg.det(np.block([[r2-r0], [r12]]))
    # Determine if all conditions are met, allowing for a small negative
    # tolerance
    inside=all(cond > -tol)
    # Calculate residual of conditions
    res=np.min(cond)

    return inside,res

def find_intersections_mat(Xp1,Yp1,Zp1,Xp2,Yp2,Zp2):
    # Calculates intersections between the spanwise segments of surface
    # It is a Python alternative to the C function find_intersectionsso.c but
    # it is much slower.
    # Xp1,Yp1,Zp1 and the panels of surface Xp2,Yp2,Zp2.
    # Xp1,Yp1,Zp1: m*n array of wing or wake vertices
    # Xp2,Yp2,Zp2: m*n array of fuselage vertices
    # xyzint: x*3 array of intersections. If surface Xp1,Yp1,Zp1 is a wing
    #     x=m+1. If it is a wake x can take any value, depending on the
    #     geometries.
    # nint: number of intersections. This can be higher than x.

    m1,n1=Xp1.shape
    m1=m1-1
    n1=n1-1
    m2,n2=Xp2.shape
    m2=m2-1
    n2=n2-1
    
    # Initialize array to hold coordinates of intersection points
    xyzint=100000.0*np.ones([m1+1,3])
    # Initialize number of intersections
    nint=0
    # Cycle through the spanwise edges of body 1
    for i1 in range (0,m1+1):
        for j1 in range (0,n1):
            # Line segment of body 1
            x0=Xp1[i1,j1:j1+2]
            y0=Yp1[i1,j1:j1+2]
            z0=Zp1[i1,j1:j1+2]
            # Cycle through all panels of body 2
            for i2 in range (0,m2):
                for j2 in range (0,n2):
                    # Split each panel into two planar triangles
                    # Start with first triangle
                    x=np.array([Xp2[i2,j2], Xp2[i2,j2+1], Xp2[i2+1,j2+1]])
                    y=np.array([Yp2[i2,j2], Yp2[i2,j2+1], Yp2[i2+1,j2+1]])
                    z=np.array([Zp2[i2,j2], Zp2[i2,j2+1], Zp2[i2+1,j2+1]])
                    # Calculate area of triangle
                    area=1.0/2.0*(x[0]*(y[1]-y[2])+x[1]*(y[2]-y[0])+x[2]*(y[0]-y[1]))
                    # The nose and tail of the fuselage only have one triangle
                    # but four vertices. This means that one of the two
                    # triangles has zero area. Exclude this from the
                    # calculation.
                    if area != 0:        
                        # Calculate control point of current triangle
                        xc=np.mean(x)
                        yc=np.mean(y)
                        zc=np.mean(z)
                        # Calculate unit vector normal to the triangle's plane
                        normal=np.cross(np.array([x[2]-x[0], y[2]-y[0], z[2]-z[0]]),np.array([x[1]-x[0], y[1]-y[0], z[1]-z[0]]))
                        normal=normal/np.sqrt(np.sum(normal*normal))
                        # Calculate first tangent unit vector parallel to side 1
                        tau1=np.array([x[1]-x[0], y[1]-y[0], z[1]-z[0]])
                        tau1=tau1/np.sqrt(np.sum(tau1*tau1))
                        # Calculate second tangent unit vector orthogonal to normal and tau1
                        tau2=np.cross(normal,tau1)  
                        # Calculate coordinates of panel vertices in local panel coordinates
                        xl=np.zeros(3)
                        yl=np.zeros(3)
                        zl=np.zeros(3)
                        for iv in range (0,3):
                            xl[iv]=(x[iv]-xc)*tau1[0]+(y[iv]-yc)*tau1[1]+(z[iv]-zc)*tau1[2]
                            yl[iv]=(x[iv]-xc)*tau2[0]+(y[iv]-yc)*tau2[1]+(z[iv]-zc)*tau2[2]
                            zl[iv]=(x[iv]-xc)*normal[0]+(y[iv]-yc)*normal[1]+(z[iv]-zc)*normal[2]
                        # End for
                        # Calculate coordinates of line segment in panel's local frame of
                        # reference
                        x0l=(x0-xc)*tau1[0]+(y0-yc)*tau1[1]+(z0-zc)*tau1[2]
                        y0l=(x0-xc)*tau2[0]+(y0-yc)*tau2[1]+(z0-zc)*tau2[2]
                        z0l=(x0-xc)*normal[0]+(y0-yc)*normal[1]+(z0-zc)*normal[2]
                        # Calculate unit vector tangent to line passing through wing points
                        tauw=np.array([x0l[1]-x0l[0], y0l[1]-y0l[0], z0l[1]-z0l[0]])
                        tauw=tauw/np.sqrt(np.sum(tauw*tauw))
                        # Calculate intersection between extension of line
                        # x0l,y0l,z0l and plane at zl=0.
                        # Equation of line is [x0l(1) y0l(1) z0l(1)]+sat0*tauw=[xat0 yat0 0]
                        # Third equation: z0l(1)+sat0*tauw(3)=0
                        sat0=-z0l[0]/tauw[2]
                        # First equation
                        xat0=x0l[0]+sat0*tauw[0]
                        # Second equation
                        yat0=y0l[0]+sat0*tauw[1]
                        # Determine if point (xat0,yat0) lies inside current panel
                        inside,_=in_triangle(xl,yl,xat0,yat0,1e-12)
                        if inside and np.prod(z0l) < 0:     # Only counts as an intersection if the signs of z0l are different
                            nint=nint+1
                            # Calculate coordinates of intersection point in global
                            # coordinates
                            xat0g=tau1[0]*xat0+tau2[0]*yat0+xc
                            yat0g=tau1[1]*xat0+tau2[1]*yat0+yc
                            zat0g=tau1[2]*xat0+tau2[2]*yat0+zc
                            # Store coordinates of intersection point
                            xyzint[i1,:]=np.array([xat0g, yat0g, zat0g]) 
                        # End if inside
                    # End if area
                    
                    if not inside:
                        # If we have not found an intersection yet, check the second triangle
                        x=np.array([Xp2[i2+1,j2+1], Xp2[i2+1,j2], Xp2[i2,j2]])
                        y=np.array([Yp2[i2+1,j2+1], Yp2[i2+1,j2], Yp2[i2,j2]])
                        z=np.array([Zp2[i2+1,j2+1], Zp2[i2+1,j2], Zp2[i2,j2]])
                        # Calculate area of triangle
                        area=1.0/2.0*(x[0]*(y[1]-y[2])+x[1]*(y[2]-y[0])+x[2]*(y[0]-y[1]))
                        # The nose and tail of the fuselage only have one triangle
                        # but four vertices. This means that one of the two
                        # triangles has zero area. Exclude this from the
                        # calculation.
                        if area != 0:        
                            # Calculate control point of current triangle
                            xc=np.mean(x)
                            yc=np.mean(y)
                            zc=np.mean(z)
                            # Calculate unit vector normal to the triangle's plane
                            normal=np.cross(np.array([x[2]-x[0], y[2]-y[0], z[2]-z[0]]),np.array([x[1]-x[0], y[1]-y[0], z[1]-z[0]]))
                            normal=normal/np.sqrt(np.sum(normal*normal))
                            # Calculate first tangent unit vector parallel to side 1
                            tau1=np.array([x[1]-x[0], y[1]-y[0], z[1]-z[0]])
                            tau1=tau1/np.sqrt(np.sum(tau1*tau1))
                            # Calculate second tangent unit vector orthogonal to normal and tau1
                            tau2=np.cross(normal,tau1)  
                            # Calculate coordinates of panel vertices in local panel coordinates
                            xl=np.zeros(3)
                            yl=np.zeros(3)
                            zl=np.zeros(3)
                            for iv in range (0,3):
                                xl[iv]=(x[iv]-xc)*tau1[0]+(y[iv]-yc)*tau1[1]+(z[iv]-zc)*tau1[2]
                                yl[iv]=(x[iv]-xc)*tau2[0]+(y[iv]-yc)*tau2[1]+(z[iv]-zc)*tau2[2]
                                zl[iv]=(x[iv]-xc)*normal[0]+(y[iv]-yc)*normal[1]+(z[iv]-zc)*normal[2]
                            # End for
                            # Calculate coordinates of line segment in panel's local frame of
                            # reference
                            x0l=(x0-xc)*tau1[0]+(y0-yc)*tau1[1]+(z0-zc)*tau1[2]
                            y0l=(x0-xc)*tau2[0]+(y0-yc)*tau2[1]+(z0-zc)*tau2[2]
                            z0l=(x0-xc)*normal[0]+(y0-yc)*normal[1]+(z0-zc)*normal[2]
                            # Calculate unit vector tangent to line passing through wing points
                            tauw=np.array([x0l[1]-x0l[0], y0l[1]-y0l[0], z0l[1]-z0l[0]])
                            tauw=tauw/np.sqrt(np.sum(tauw*tauw))
                            # Calculate intersection between extension of line
                            # x0l,y0l,z0l and plane at zl=0.
                            # Equation of line is [x0l(1) y0l(1) z0l(1)]+sat0*tauw=[xat0 yat0 0]
                            # Third equation: z0l(1)+sat0*tauw(3)=0
                            sat0=-z0l[0]/tauw[2]
                            # First equation
                            xat0=x0l[0]+sat0*tauw[0]
                            # Second equation
                            yat0=y0l[0]+sat0*tauw[1]
                            # Determine if point (xat0,yat0) lies inside current panel
                            inside,_=in_triangle(xl,yl,xat0,yat0,1e-12)
                            if inside and np.prod(z0l) < 0:     # Only counts as an intersection if the signs of z0l are different
                                nint=nint+1
                                # Calculate coordinates of intersection point in global
                                # coordinates
                                xat0g=tau1[0]*xat0+tau2[0]*yat0+xc
                                yat0g=tau1[1]*xat0+tau2[1]*yat0+yc
                                zat0g=tau1[2]*xat0+tau2[2]*yat0+zc
                                # Store coordinates of intersection point
                                xyzint[i1,:]=np.array([xat0g, yat0g, zat0g]) 
                            # End if inside
                        # End if area  
                    # End if inside
                # End for j2
            # End for i2
        # End for j1
    # End for i1
 
    # Truncate xyzint if the numnber of intersections is lower than the first
    # dimension of xyzint
    if nint < xyzint.shape[0]:
        xyzint=xyzint[0:nint,:]
    # Enf if

    return xyzint,nint

def Cfuncs_find_intersections(install_dir):
    # Loads the c function used for finding intersections between wings and bodies. 
    # install_dir: Path to ./Common directory
    # find_intersections: function to calculate intersections between wings and bodies
    
    # Load c function 
    if platform.system() == 'Windows':
        so_file = install_dir+r"\find_intersectionsso.so" 
    else:
        so_file = install_dir+r"find_intersectionsso.so" 
    # End if
    lib = ctypes.cdll.LoadLibrary(so_file)
    find_intersections = lib.find_intersections
    find_intersections.restype = None
    # Define arguments of sdpminf
    find_intersections.argtypes = [ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
    		        ctypes.c_int,ctypes.c_int,ctypes.c_int,ctypes.c_int,
                    ctypes.c_int,ctypes.c_int]
    
    return find_intersections
        
def wing_fuselage_intersections(body,ibody1,ibody2,j1min,install_dir):
    # Find intersections between a horizontal wing surface and a fuselage.
    # Truncate the wing to the intersections and split the fuselage to lower
    # and upper sections, i.e. below the wing and above the wing.
    # body: struct array containing the geometric information of all bodies.
    #       The length of body must be at least 2 (a wing and a fuselage).
    #       NB: struct body has one additional element after execution.
    #       body(ibody1) remains the wing. body(ibody2) becomes the lower
    #       fuselage and body(ibody2+1) becomes the upper fuselage. Ensure that
    #       body did not contain an ibody2+1 element before execution because
    #       it will be overwritten.
    # ibody1: Index of wing in body array.
    # ibody2: Index of fuselage in body array.
    # jmin: Spanwise index of wing from which to start looking for
    #       intersections.
    # xyzint: (body(ibody1).m+1)*3 array of coordinates of intersections.
    
    # Load c function to calculate the steady SDPM influence coefficient
    # matrices
    find_intersections=Cfuncs_find_intersections(install_dir)
    
    # Find intersections between wing surface in ibody1 and fuselage surface in
    # ibody2
    # Calculate maximum spanwise index of body 1 that must be included in the
    # search for intersections.
    j1max=np.argwhere(body['Yp0'][ibody1][0,:] > np.max(body['Yp0'][ibody2]))
    j1max=j1max[0][0]
    # Extract x, y and z coordinates of panel vertices between j1min, j1max and
    # j2min, j2max
    Xp1=np.ascontiguousarray(body['Xp0'][ibody1][:,j1min:j1max+1])
    Yp1=np.ascontiguousarray(body['Yp0'][ibody1][:,j1min:j1max+1])
    Zp1=np.ascontiguousarray(body['Zp0'][ibody1][:,j1min:j1max+1])
    Xp2=np.ascontiguousarray(body['Xp0'][ibody2])
    Yp2=np.ascontiguousarray(body['Yp0'][ibody2])
    Zp2=np.ascontiguousarray(body['Zp0'][ibody2])
    xyzint1=100000.0*np.ones([Xp1.shape[0],3])
    # Calculate intersections
    find_intersections(xyzint1,Xp1,Yp1,Zp1,Xp2,Yp2,Zp2,Xp1.shape[0],Xp1.shape[1],Xp2.shape[0],Xp2.shape[1],xyzint1.shape[0],xyzint1.shape[1])
    # Determine number of intersections
    iko=np.argwhere(xyzint1[:,0] != 100000.0)
    nint1=len(iko)
    xyzint1=xyzint1[iko[:,0],:]    
    
    # Find intersections between wake surface in ibody1 and fuselage surface in
    # ibody2
    # j1max is expected to be the same as before because the spanwise wake
    # positions coincide with those of the wing's surface.
    # Calculate intersections
    Xp1=np.ascontiguousarray(body['Xw0'][ibody1][:,j1min:j1max+1])
    Yp1=np.ascontiguousarray(body['Yw0'][ibody1][:,j1min:j1max+1])
    Zp1=np.ascontiguousarray(body['Zw0'][ibody1][:,j1min:j1max+1])
    xyzint2=100000.0*np.ones([Xp1.shape[0],3])
    find_intersections(xyzint2,Xp1,Yp1,Zp1,Xp2,Yp2,Zp2,Xp1.shape[0],Xp1.shape[1],Xp2.shape[0],Xp2.shape[1],xyzint2.shape[0],xyzint2.shape[1])
    # Determine number of intersections
    iko=np.argwhere(xyzint2[:,0] != 100000.0)
    nint2=len(iko)
    xyzint2=xyzint2[iko[:,0],:] 
    
    # The tail of the fuselage will, in general, not lie on an interception with
    # the wake. Calculate the interecption of the tail of the fuselage with the
    # wake.
    # Isolate the last column of the fuselage (tail section)
    Xp1=np.ascontiguousarray(np.reshape(body['Xp0'][ibody2][:,body['n'][ibody2]],(1,body['m'][ibody2]+1)))
    Yp1=np.ascontiguousarray(np.reshape(body['Yp0'][ibody2][:,body['n'][ibody2]],(1,body['m'][ibody2]+1)))
    Zp1=np.ascontiguousarray(np.reshape(body['Zp0'][ibody2][:,body['n'][ibody2]],(1,body['m'][ibody2]+1)))
    Xp2=np.ascontiguousarray(body['Xw0'][ibody1])
    Yp2=np.ascontiguousarray(body['Yw0'][ibody1])
    Zp2=np.ascontiguousarray(body['Zw0'][ibody1])
    xyzint3=100000.0*np.ones([Xp1.shape[0],3])
    find_intersections(xyzint3,Xp1,Yp1,Zp1,Xp2,Yp2,Zp2,Xp1.shape[0],Xp1.shape[1],Xp2.shape[0],Xp2.shape[1],xyzint3.shape[0],xyzint3.shape[1])
    # Determine number of intersections
    iko=np.argwhere(xyzint3[:,0] != 100000.0)
    nint3=len(iko)
    xyzint3=xyzint3[iko[:,0],:] 
    
    # Truncate the wing to start at the fuselage intersections
    Xp0=np.block([np.reshape(xyzint1[:,0],(body['m'][ibody1]+1,1)), body['Xp0'][ibody1][:,j1max:body['n'][ibody1]+1]])
    Yp0=np.block([np.reshape(xyzint1[:,1],(body['m'][ibody1]+1,1)), body['Yp0'][ibody1][:,j1max:body['n'][ibody1]+1]])
    Zp0=np.block([np.reshape(xyzint1[:,2],(body['m'][ibody1]+1,1)), body['Zp0'][ibody1][:,j1max:body['n'][ibody1]+1]])

    if nint2 > 1 and nint3 > 0: # If there are interceptions between the wake and the fuselage all the way from the trailing edge to the tail
        # Find index of wake outboard from last intersection with fuselage
        iy=np.argwhere(body['Yw0'][ibody1][0,:] > xyzint3[0,1])
        iy=iy[0][0]
        # Find index of wake at iy downstream from last intersection with fuselage
        ix=np.argwhere(body['Xw0'][ibody1][:,iy] > xyzint3[0,0])
        ix=ix[0][0]
        # Move x-coordinate of wake index ix to obtain row at last intersection with fuselage
        dx=xyzint3[0,0]-xyzint2[ix-1,0] #  Distance between last two interceptions between the wake and the fusalge
        body['Xw0'][ibody1][ix,:]=body['Xw0'][ibody1][ix-1,:]+dx
        
        # Truncate the wake to start at the fuselage intersections
        nint=ix-1   # Number of intersections
        Xw0=body['Xw0'][ibody1][:,j1max-1:body['n'][ibody1]+1]  # Truncate x coordinates of wake
        Xw0[0:nint+1,0]=xyzint2[0:nint+1,0]                             # Assign intersections with fuselage at first row of wake
        Xw0[nint+1,0]=xyzint3[0,0]                                      # Assign fuselage tail at first row of wake
        Xw0[nint+2:body['mw'][ibody1]+1,0]=xyzint3[0,0]+Xw0[nint+2:body['mw'][ibody1]+1,1]-Xw0[nint+1,1] # Increment x coordinate of first fow using increments from second row
        Yw0=body['Yw0'][ibody1][:,j1max-1:body['n'][ibody1]+1]  # Truncate y coordinates of wake
        Yw0[0:nint+1,0]=xyzint2[0:nint+1,1]                             # Assign intersections with fuselage at first row of wake
        Yw0[nint+1:body['mw'][ibody1]+1,0]=xyzint3[0,1]             # Assign fuselage tail at first row of wake
        Zw0=body['Zw0'][ibody1][:,j1max-1:body['n'][ibody1]+1]  # Truncate y coordinates of wake
        Zw0[0:nint+1,0]=xyzint2[0:nint+1,2]                             # Assign intersections with fuselage at first row of wake
        Zw0[nint+1:body['mw'][ibody1]+1,0]=xyzint3[0,2]             # Assign fuselage tail at first row of wake
        # Calculate new x increments
        wake_dx=Xw0[1:body['mw'][ibody1]+1,0]-Xw0[0:body['mw'][ibody1],0]
    else:   # If the only interception between the wake and the fuselage lies at the wing root
        dx=np.reshape(np.block([0.0,np.cumsum(body['wake_dx'][ibody1])]),(body['mw'][ibody1]+1,1))
        Xw0=np.block([xyzint1[0,0]+dx,body['Xw0'][ibody1][:,j1max:body['n'][ibody1]+1]])
        Yw0=np.block([xyzint1[0,1]*np.ones([body['mw'][ibody1]+1,1]),body['Yw0'][ibody1][:,j1max:body['n'][ibody1]+1]])
        Zw0=np.block([xyzint1[0,2]*np.ones([body['mw'][ibody1]+1,1]),body['Zw0'][ibody1][:,j1max:body['n'][ibody1]+1]])
        wake_dx=body['wake_dx'][ibody1]
    # End if

    # Assign wing
    m,n=Xp0.shape
    m=m-1
    n=n-1

    # Calculate control points, normal vectors and areas
    Xc0,Yc0,Zc0,nx0,ny0,nz0,tauxx0,tauxy0,tauxz0,tauyx0,tauyy0,tauyz0,s0,cpln0,sm0,tmx0,tmy0,tmz0,sn0,tnx0,tny0,tnz0=SDPM_control_normal_tangential(Xp0,Yp0,Zp0,body['dir_tau'][ibody1])

    # Reshape control point coordinate matrices into vectors
    Xc0all=np.reshape(Xc0,(m*n,1),order='C')
    Yc0all=np.reshape(Yc0,(m*n,1),order='C')
    Zc0all=np.reshape(Zc0,(m*n,1),order='C')
    # Reshape normal vector component matrices into vectors
    nx0all=np.reshape(nx0,(m*n,1),order='C')
    ny0all=np.reshape(ny0,(m*n,1),order='C')
    nz0all=np.reshape(nz0,(m*n,1),order='C')
    # Reshape panel area matrix into a vector
    s0all=np.reshape(s0,(m*n,1),order='C')

    # Assign data to body struct array
    body['Xp0'][ibody1]=Xp0
    body['Yp0'][ibody1]=Yp0
    body['Zp0'][ibody1]=Zp0
    body['Xc0'][ibody1]=Xc0
    body['Yc0'][ibody1]=Yc0
    body['Zc0'][ibody1]=Zc0 
    body['Xc0all'][ibody1]=Xc0all
    body['Yc0all'][ibody1]=Yc0all
    body['Zc0all'][ibody1]=Zc0all 
    body['Xw0'][ibody1]=Xw0
    body['Yw0'][ibody1]=Yw0
    body['Zw0'][ibody1]=Zw0
    body['mw'][ibody1]=body['mw'][ibody1]
    body['wake_dx'][ibody1]=wake_dx
    body['nx0'][ibody1]=nx0
    body['ny0'][ibody1]=ny0
    body['nz0'][ibody1]=nz0
    body['tauxx0'][ibody1]=tauxx0
    body['tauxy0'][ibody1]=tauxy0
    body['tauxz0'][ibody1]=tauxz0
    body['tauyx0'][ibody1]=tauyx0
    body['tauyy0'][ibody1]=tauyy0
    body['tauyz0'][ibody1]=tauyz0
    body['nx0all'][ibody1]=nx0all
    body['ny0all'][ibody1]=ny0all
    body['nz0all'][ibody1]=nz0all
    body['s0'][ibody1]=s0
    body['s0all'][ibody1]=s0all
    body['m'][ibody1]=m
    body['n'][ibody1]=n
    body['c0'][ibody1]=body['c0'][ibody1]
    body['b'][ibody1]=body['b'][ibody1]
    body['yc'][ibody1]=body['yc'][ibody1]
    body['AR'][ibody1]=body['AR'][ibody1]
    body['S'][ibody1]=body['S'][ibody1]
    body['name'][ibody1]=body['name'][ibody1]
    body['dir_tau'][ibody1]=body['dir_tau'][ibody1]
    
    xyzint=xyzint1
    # Indices of fuselage sections aft of the wing trailing edge
    ite=np.argwhere(body['Xp0'][ibody2][0,:] > xyzint[0,0])
    # Indices of fuselage sections forward of the wing leading edge
    ile=np.argwhere(body['Xp0'][ibody2][0,:] < xyzint[body['m'][ibody1]//2,0])

    # Create lower fuselage
    # Determine number of heightwise vertices on lower fuselage by looking at
    # how many vertices lie lower than the leading edge intersection.
    ilo=np.argwhere(body['Zp0'][ibody2][:,ile[-1,0]+1] < xyzint[body['m'][ibody1]//2,2])
    ilo=np.max(ilo)+1   # Add one for the intersection point
    # Resample lower fuselage at the wing root
    Xp_int=np.zeros([ilo+1,body['m'][ibody1]//2+1])
    Yp_int=np.zeros([ilo+1,body['m'][ibody1]//2+1])
    Zp_int=np.zeros([ilo+1,body['m'][ibody1]//2+1])
    theta_res=np.zeros(body['m'][ibody1]//2+1)
    for i in range (body['m'][ibody1]//2,-1,-1):
        Xp_int[:,body['m'][ibody1]//2-i]=xyzint[i,0]*np.ones(ilo+1)
        iko=np.argwhere(body['Xp0'][ibody2][0,:] < xyzint[i,0])
        jhere=np.array([np.max(iko), np.max(iko)+1])
        xc=body['Xp0'][ibody2][0,jhere]
        # Interpolate linearly between sections
        slope=(body['Yp0'][ibody2][:,jhere[1]]-body['Yp0'][ibody2][:,jhere[0]])/(xc[1]-xc[0])
        Y=slope*(xyzint[i,0]-xc[0])+body['Yp0'][ibody2][:,jhere[0]]
        slope=(body['Zp0'][ibody2][:,jhere[1]]-body['Zp0'][ibody2][:,jhere[0]])/(xc[1]-xc[0])
        Z=slope*(xyzint[i,0]-xc[0])+body['Zp0'][ibody2][:,jhere[0]]
        # Calculate midpoint of z coordinates
        Zmid=(Z[-1]-Z[0])/2.0+Z[0]
        # Angle of intersection point on lower wing
        theta_res[body['m'][ibody1]//2-i]=np.arctan2(xyzint[i,2]-Zmid,xyzint[i,1])
        # Transform Y and Z to polar coordinates
        th=np.arctan2(Z-Zmid,Y)
        r=np.sqrt(Y**2+(Z-Zmid)**2)   
        # All angles from bottom of fuselage to intersection point
        theta=np.linspace(np.min(th),theta_res[body['m'][ibody1]//2-i],ilo+1)
        # Interpolate onto angles theta
        cs = CubicSpline(th, r)
        # Calculate airfoil shape at xp
        rint=cs(theta)
        # Transform to Cartesian coordinates
        Yint=rint*np.cos(theta)
        Zint=rint*np.sin(theta)
        # Store resampled y and z coordinates
        Yp_int[:,body['m'][ibody1]//2-i]=Yint
        Zp_int[:,body['m'][ibody1]//2-i]=Zint+Zmid
        # Ensure that last point is the intersection point
        Yp_int[ilo,body['m'][ibody1]//2-i]=xyzint[i,1]
        Zp_int[ilo,body['m'][ibody1]//2-i]=xyzint[i,2]
    # End for
    
    # Resample lower fuselage forward of the wing leading edge
    Xp_for=np.zeros([ilo+1,len(ile)])
    Yp_for=np.zeros([ilo+1,len(ile)])
    Zp_for=np.zeros([ilo+1,len(ile)])
    # Set coordinates of nose section
    Xp_for[:,0]=body['Xp0'][ibody2][0,0]
    Zp_for[:,0]=body['Zp0'][ibody2][0,0]
    # Cycle for all other sections
    for i in range (1,len(ile)):
        Xp_for[:,i]=body['Xp0'][ibody2][0,ile[i]]
        Y=body['Yp0'][ibody2][:,ile[i]]
        Z=body['Zp0'][ibody2][:,ile[i]]
        # Calculate midpoint of z coordinates
        Zmid=(Z[-1,0]-Z[0,0])/2.0+Z[0,0]
        # Transform Y and Z to polar coordinates
        th=np.arctan2(Z[:,0]-Zmid,Y[:,0])
        r=np.sqrt(Y[:,0]**2+(Z[:,0]-Zmid)**2)   
        # All angles from bottom of fuselage to intersection point
        theta=np.linspace(np.min(th),theta_res[0],ilo+1)
        # Interpolate onto angles theta
        cs = CubicSpline(th, r)
        # Calculate airfoil shape at xp
        rint=cs(theta)
        # Transform to Cartesian coordinates
        Yint=rint*np.cos(theta)
        Zint=rint*np.sin(theta)
        # Store resampled y and z coordinates
        Yp_for[:,i]=Yint
        Zp_for[:,i]=Zint+Zmid       
    # End for

    # Resample lower fuselage aft of the wing trailing edge
    Xp_aft=np.zeros([ilo+1,len(ite)])
    Yp_aft=np.zeros([ilo+1,len(ite)])
    Zp_aft=np.zeros([ilo+1,len(ite)])
    # Cycle for all sections
    for i in range (0,len(ite)):
        Xp_aft[:,i]=body['Xp0'][ibody2][0,ite[i]]
        Y=body['Yp0'][ibody2][:,ite[i]]
        Z=body['Zp0'][ibody2][:,ite[i]]
        # Calculate midpoint of z coordinates
        Zmid=(Z[-1,0]-Z[0,0])/2.0+Z[0,0]
        # Transform Y and Z to polar coordinates
        th=np.arctan2(Z[:,0]-Zmid,Y[:,0])
        r=np.sqrt(Y[:,0]**2+(Z[:,0]-Zmid)**2)   
        # All angles from bottom of fuselage to intersection point
        theta=np.linspace(np.min(th),theta_res[0],ilo+1)
        # Interpolate onto angles theta
        cs = CubicSpline(th, r)
        # Calculate airfoil shape at xp
        rint=cs(theta)
        # Transform to Cartesian coordinates
        Yint=rint*np.cos(theta)
        Zint=rint*np.sin(theta)
        # Store resampled y and z coordinates
        Yp_aft[:,i]=Yint
        Zp_aft[:,i]=Zint+Zmid       
    # End for

    if nint2 > 1 and nint3 > 0: # If there are interceptions between the wake and the fuselage all the way from the trailing edge to the tail
        # Now interpolate lower aft fuselage to the wake intercepts
        xyzint2=np.block([[xyzint2[1:nint+1,:]],[xyzint3]])     # Assemble all wake intercepts
        # Initialize interpolated aft fuselage arrays
        Xp_aft_int=np.zeros([ilo+1,nint+1])
        Yp_aft_int=np.zeros([ilo+1,nint+1])
        Zp_aft_int=np.zeros([ilo+1,nint+1])
        # Cycle through the lengthwise sections
        for i in range (0,nint+1):
            # Assign value of x at this lengthwise section
            Xp_aft_int[:,i]=xyzint2[i,0]*np.ones(ilo+1)
            # Find which index of Xp_aft is higher than the current x intercept.
            iko=np.argwhere(Xp_aft[0,:] < xyzint2[i,0])
            if len(iko) != 0:
                jhere=np.array([np.max(iko), np.max(iko)+1])
                xc=Xp_aft[0,jhere]
                # Interpolate linearly between sections
                slope=(Yp_aft[:,jhere[1]]-Yp_aft[:,jhere[0]])/(xc[1]-xc[0])
                Y=slope*(xyzint2[i,0]-xc[0])+Yp_aft[:,jhere[0]]
                slope=(Zp_aft[:,jhere[1]]-Zp_aft[:,jhere[0]])/(xc[1]-xc[0])
                Z=slope*(xyzint2[i,0]-xc[0])+Zp_aft[:,jhere[0]]
            else:   # If iko is empty, we need to use the last section of Xp_int, Yp_int, Zp_int
                xc=np.array([Xp_int[0,-1], Xp_aft[0,0]])
                slope=(Yp_aft[:,0]-Yp_int[:,-1])/(xc[1]-xc[0])
                Y=slope*(xyzint2[i,0]-xc[0])+Yp_int[:,-1]
                slope=(Zp_aft[:,0]-Zp_int[:,-1])/(xc[1]-xc[0])
                Z=slope*(xyzint2[i,0]-xc[0])+Zp_int[:,-1]
            # End if
            Yp_aft_int[:,i]=Y
            Zp_aft_int[:,i]=Z
        # End for
        # Ensure that top of fuselage lies on wake intercepts
        Yp_aft_int[ilo,:]=xyzint2[:,1]
        Zp_aft_int[ilo,:]=xyzint2[:,2]
        # Assemble coordinates of lower fuselage
        Xp_lo=np.block([Xp_for, Xp_int, Xp_aft_int])
        Yp_lo=np.block([Yp_for, Yp_int, Yp_aft_int])
        Zp_lo=np.block([Zp_for, Zp_int, Zp_aft_int])
    else:   # If the only interception between the wake and the fuselage lies at the wing root
        # Assemble coordinates of lower fuselage
        Xp_lo=np.block([Xp_for, Xp_int, Xp_aft])
        Yp_lo=np.block([Yp_for, Yp_int, Yp_aft])
        Zp_lo=np.block([Zp_for, Zp_int, Zp_aft])
    # End if
    
    # Create upper fuselage
    # Set number of heightwise vertices on upper fuselage
    iup=body['m'][ibody2]-(ilo-1)
    
    # Resample upper fuselage at the wing root
    Xp_int=np.zeros([iup+1,body['m'][ibody1]//2+1])
    Yp_int=np.zeros([iup+1,body['m'][ibody1]//2+1])
    Zp_int=np.zeros([iup+1,body['m'][ibody1]//2+1])
    for i in range (0,body['m'][ibody1]//2+1):
        Xp_int[:,i]=xyzint[i+body['m'][ibody1]//2,0]*np.ones(iup+1)
        iko=np.argwhere(body['Xp0'][ibody2][0,:] < xyzint[i+body['m'][ibody1]//2,0])
        jhere=np.array([np.max(iko), np.max(iko)+1])
        xc=body['Xp0'][ibody2][0,jhere]
        # Interpolate linearly between sections
        slope=(body['Yp0'][ibody2][:,jhere[1]]-body['Yp0'][ibody2][:,jhere[0]])/(xc[1]-xc[0])
        Y=slope*(xyzint[i+body['m'][ibody1]//2,0]-xc[0])+body['Yp0'][ibody2][:,jhere[0]]
        slope=(body['Zp0'][ibody2][:,jhere[1]]-body['Zp0'][ibody2][:,jhere[0]])/(xc[1]-xc[0])
        Z=slope*(xyzint[i+body['m'][ibody1]//2,0]-xc[0])+body['Zp0'][ibody2][:,jhere[0]]
        # Calculate midpoint of z coordinates
        Zmid=(Z[-1]-Z[0])/2.0+Z[0]
        # Angle of intersection point on lower wing
        theta_res[i]=np.arctan2(xyzint[i+body['m'][ibody1]//2,2]-Zmid,xyzint[i+body['m'][ibody1]//2,1])
        # Transform Y and Z to polar coordinates
        th=np.arctan2(Z-Zmid,Y)
        r=np.sqrt(Y**2+(Z-Zmid)**2)   
        # All angles from bottom of fuselage to intersection point
        theta=np.linspace(theta_res[i],np.max(th),iup+1)
        # Interpolate onto angles theta
        cs = CubicSpline(th, r)
        # Calculate airfoil shape at xp
        rint=cs(theta)
        # Transform to Cartesian coordinates
        Yint=rint*np.cos(theta)
        Zint=rint*np.sin(theta)
        # Store resampled y and z coordinates
        Yp_int[:,i]=Yint
        Zp_int[:,i]=Zint+Zmid
        # Ensure that last point is the intersection point
        Yp_int[0,i]=xyzint[i+body['m'][ibody1]//2,1]
        Zp_int[0,i]=xyzint[i+body['m'][ibody1]//2,2]
    # End for    

    # Resample upper fuselage forward of the wing leading edge
    Xp_for=np.zeros([iup+1,len(ile)])
    Yp_for=np.zeros([iup+1,len(ile)])
    Zp_for=np.zeros([iup+1,len(ile)])
    # Set coordinates of nose section
    Xp_for[:,0]=body['Xp0'][ibody2][0,0]
    Zp_for[:,0]=body['Zp0'][ibody2][0,0]
    # Cycle for all other sections
    for i in range (1,len(ile)):
        Xp_for[:,i]=body['Xp0'][ibody2][0,ile[i]]
        Y=body['Yp0'][ibody2][:,ile[i]]
        Z=body['Zp0'][ibody2][:,ile[i]]
        # Calculate midpoint of z coordinates
        Zmid=(Z[-1,0]-Z[0,0])/2.0+Z[0,0]
        # Transform Y and Z to polar coordinates
        th=np.arctan2(Z[:,0]-Zmid,Y[:,0])
        r=np.sqrt(Y[:,0]**2+(Z[:,0]-Zmid)**2)   
        # All angles from bottom of fuselage to intersection point
        theta=np.linspace(theta_res[0],np.max(th),iup+1)
        # Interpolate onto angles theta
        cs = CubicSpline(th, r)
        # Calculate airfoil shape at xp
        rint=cs(theta)
        # Transform to Cartesian coordinates
        Yint=rint*np.cos(theta)
        Zint=rint*np.sin(theta)
        # Store resampled y and z coordinates
        Yp_for[:,i]=Yint
        Zp_for[:,i]=Zint+Zmid       
    # End for

    # Resample upper fuselage aft of the wing trailing edge
    Xp_aft=np.zeros([iup+1,len(ite)])
    Yp_aft=np.zeros([iup+1,len(ite)])
    Zp_aft=np.zeros([iup+1,len(ite)])
    # Cycle for all sections
    for i in range (0,len(ite)):
        Xp_aft[:,i]=body['Xp0'][ibody2][0,ite[i]]
        Y=body['Yp0'][ibody2][:,ite[i]]
        Z=body['Zp0'][ibody2][:,ite[i]]
        # Calculate midpoint of z coordinates
        Zmid=(Z[-1,0]-Z[0,0])/2.0+Z[0,0]
        # Transform Y and Z to polar coordinates
        th=np.arctan2(Z[:,0]-Zmid,Y[:,0])
        r=np.sqrt(Y[:,0]**2+(Z[:,0]-Zmid)**2)   
        # All angles from bottom of fuselage to intersection point
        theta=np.linspace(theta_res[body['m'][ibody1]//2],np.max(th),iup+1)
        # Interpolate onto angles theta
        cs = CubicSpline(th, r)
        # Calculate airfoil shape at xp
        rint=cs(theta)
        # Transform to Cartesian coordinates
        Yint=rint*np.cos(theta)
        Zint=rint*np.sin(theta)
        # Store resampled y and z coordinates
        Yp_aft[:,i]=Yint
        Zp_aft[:,i]=Zint+Zmid       
    # End for
        
    if nint2 > 1 and nint3 > 0: # If there are interceptions between the wake and the fuselage all the way from the trailing edge to the tail
        # Now interpolate upper aft fuselage to the wake intercepts
        # Initialize interpolated aft fuselage arrays
        Xp_aft_int=np.zeros([iup+1,nint+1])
        Yp_aft_int=np.zeros([iup+1,nint+1])
        Zp_aft_int=np.zeros([iup+1,nint+1])
        # Cycle through the lengthwise sections
        for i in range (0,nint+1):
            # Assign value of x at this lengthwise section
            Xp_aft_int[:,i]=xyzint2[i,0]*np.ones(iup+1)
            # Find which index of Xp_aft is higher than the current x intercept.
            iko=np.argwhere(Xp_aft[0,:] < xyzint2[i,0])
            if len(iko) != 0:
                jhere=np.array([np.max(iko), np.max(iko)+1])
                xc=Xp_aft[0,jhere]
                # Interpolate linearly between sections
                slope=(Yp_aft[:,jhere[1]]-Yp_aft[:,jhere[0]])/(xc[1]-xc[0])
                Y=slope*(xyzint2[i,0]-xc[0])+Yp_aft[:,jhere[0]]
                slope=(Zp_aft[:,jhere[1]]-Zp_aft[:,jhere[0]])/(xc[1]-xc[0])
                Z=slope*(xyzint2[i,0]-xc[0])+Zp_aft[:,jhere[0]]
            else:   # If iko is empty, we need to use the last section of Xp_int, Yp_int, Zp_int
                xc=np.array([Xp_int[0,-1], Xp_aft[0,0]])
                slope=(Yp_aft[:,0]-Yp_int[:,-1])/(xc[1]-xc[0])
                Y=slope*(xyzint2[i,0]-xc[0])+Yp_int[:,-1]
                slope=(Zp_aft[:,0]-Zp_int[:,-1])/(xc[1]-xc[0])
                Z=slope*(xyzint2[i,0]-xc[0])+Zp_int[:,-1]
            # End if
            Yp_aft_int[:,i]=Y
            Zp_aft_int[:,i]=Z
        # End for
        # Ensure that top of fuselage lies on wake intercepts
        Yp_aft_int[0,:]=xyzint2[:,1]
        Zp_aft_int[0,:]=xyzint2[:,2]
        # Assemble coordinates of lower fuselage
        Xp_up=np.block([Xp_for, Xp_int, Xp_aft_int])
        Yp_up=np.block([Yp_for, Yp_int, Yp_aft_int])
        Zp_up=np.block([Zp_for, Zp_int, Zp_aft_int])
    else:   # If the only interception between the wake and the fuselage lies at the wing root
        # Assemble coordinates of lower fuselage
        Xp_up=np.block([Xp_for, Xp_int, Xp_aft])
        Yp_up=np.block([Yp_for, Yp_int, Yp_aft])
        Zp_up=np.block([Zp_for, Zp_int, Zp_aft])
    # End if
    
    # Ensure that forward upper and lower fuselages match up.
    Yp_up[0,ile]=Yp_lo[-1,ile]
    Zp_up[0,ile]=Zp_lo[-1,ile]
    # Ensure that aft upper and lower fuselages match up.
    Yp_up[0,ile[-1,0]+body['m'][ibody1]//2+2:Yp_up.shape[1]+1]=Yp_lo[ilo,ile[-1,0]+body['m'][ibody1]//2+2:Yp_lo.shape[1]+1]
    Zp_up[0,ile[-1,0]+body['m'][ibody1]//2+2:Zp_up.shape[1]+1]=Zp_lo[ilo,ile[-1,0]+body['m'][ibody1]//2+2:Zp_lo.shape[1]+1]
    
    # Assign lower fuselage
    m,n=Xp_lo.shape
    m=m-1
    n=n-1
    # For fuselages dir_tau=1
    dir_tau=1

    # Calculate control points, normal vectors and areas
    Xc0,Yc0,Zc0,nx0,ny0,nz0,tauxx0,tauxy0,tauxz0,tauyx0,tauyy0,tauyz0,s0,cpln0,sm0,tmx0,tmy0,tmz0,sn0,tnx0,tny0,tnz0=SDPM_control_normal_tangential(Xp_lo,Yp_lo,Zp_lo,dir_tau)

    # Reshape control point coordinate matrices into vectors
    Xc0all=np.reshape(Xc0,(m*n,1),order='C')
    Yc0all=np.reshape(Yc0,(m*n,1),order='C')
    Zc0all=np.reshape(Zc0,(m*n,1),order='C')
    # Reshape normal vector component matrices into vectors
    nx0all=np.reshape(nx0,(m*n,1),order='C')
    ny0all=np.reshape(ny0,(m*n,1),order='C')
    nz0all=np.reshape(nz0,(m*n,1),order='C')
    # Reshape panel area matrix into a vector
    s0all=np.reshape(s0,(m*n,1),order='C')

    if body['mw'][ibody2] != 0:
        # Create wake
        Xw0=np.zeros([m+1,body['mw'][ibody2]+1])       # Initialize wake doublet panel vertices
        Xw0[:,0]=Xp_lo[:,n]     # Leading edge of first row of wake doublet panels
        Xw0[:,1:body['mw'][ibody2]+1]=np.matlib.repmat(np.reshape(Xw0[:,0],(m+1,1)),1,body['mw'][ibody2])+np.matlib.repmat(np.cumsum(body['wake_dx'][ibody2]),m+1,1) # Chordwise coordinates of all other wake doublet panels
        # Set y-coordinates of all spanwise wake doublet panels to corresponding
        # spanwise y-coordinates of wing panels
        Yw0=np.matlib.repmat(np.reshape(Yp_lo[:,n],(m+1,1)),1,body['mw'][ibody2]+1)
        # Set z-coordinates of all wake  doublet panel vertices to corresponding
        # bound trailing edge spanwise vortex segment z-coordinates
        Zw0=np.matlib.repmat(np.reshape(Zp_lo[:,n],(m+1,1)),1,body['mw'][ibody2]+1)
    else:
        Xw0=0.0
        Yw0=0.0
        Zw0=0.0
    # End if
    
    # Name of body
    name_original=body['name'][ibody2].decode('UTF-8')
    name=name_original+"_lower"  

    # Assign data to body struct array
    body['Xp0'][ibody2]=Xp_lo
    body['Yp0'][ibody2]=Yp_lo
    body['Zp0'][ibody2]=Zp_lo
    body['Xc0'][ibody2]=Xc0
    body['Yc0'][ibody2]=Yc0
    body['Zc0'][ibody2]=Zc0 
    body['Xc0all'][ibody2]=Xc0all
    body['Yc0all'][ibody2]=Yc0all
    body['Zc0all'][ibody2]=Zc0all 
    body['Xw0'][ibody2]=Xw0
    body['Yw0'][ibody2]=Yw0
    body['Zw0'][ibody2]=Zw0
    body['mw'][ibody2]=body['mw'][ibody2]
    body['wake_dx'][ibody2]=body['wake_dx'][ibody2]
    body['nx0'][ibody2]=nx0
    body['ny0'][ibody2]=ny0
    body['nz0'][ibody2]=nz0
    body['tauxx0'][ibody2]=tauxx0
    body['tauxy0'][ibody2]=tauxy0
    body['tauxz0'][ibody2]=tauxz0
    body['tauyx0'][ibody2]=tauyx0
    body['tauyy0'][ibody2]=tauyy0
    body['tauyz0'][ibody2]=tauyz0
    body['nx0all'][ibody2]=nx0all
    body['ny0all'][ibody2]=ny0all
    body['nz0all'][ibody2]=nz0all
    body['s0'][ibody2]=s0
    body['s0all'][ibody2]=s0all
    body['m'][ibody2]=m
    body['n'][ibody2]=n
    body['c0'][ibody2]=body['c0'][ibody2]
    body['b'][ibody2]=body['b'][ibody2]
    body['yc'][ibody2]=body['yc'][ibody2]
    body['AR'][ibody2]=body['AR'][ibody2]
    body['S'][ibody2]=body['S'][ibody2]
    body['name'][ibody2]=name
    body['dir_tau'][ibody2]=dir_tau
    
    # Assign upper fuselage
    m,n=Xp_up.shape
    m=m-1
    n=n-1
    # For fuselages dir_tau=1
    dir_tau=1

    # Calculate control points, normal vectors and areas
    Xc0,Yc0,Zc0,nx0,ny0,nz0,tauxx0,tauxy0,tauxz0,tauyx0,tauyy0,tauyz0,s0,cpln0,sm0,tmx0,tmy0,tmz0,sn0,tnx0,tny0,tnz0=SDPM_control_normal_tangential(Xp_up,Yp_up,Zp_up,dir_tau)

    # Reshape control point coordinate matrices into vectors
    Xc0all=np.reshape(Xc0,(m*n,1),order='C')
    Yc0all=np.reshape(Yc0,(m*n,1),order='C')
    Zc0all=np.reshape(Zc0,(m*n,1),order='C')
    # Reshape normal vector component matrices into vectors
    nx0all=np.reshape(nx0,(m*n,1),order='C')
    ny0all=np.reshape(ny0,(m*n,1),order='C')
    nz0all=np.reshape(nz0,(m*n,1),order='C')
    # Reshape panel area matrix into a vector
    s0all=np.reshape(s0,(m*n,1),order='C')

    if body['mw'][ibody2] != 0:
        # Create wake
        Xw0=np.zeros([m+1,body['mw'][ibody2]+1])       # Initialize wake doublet panel vertices
        Xw0[:,0]=Xp_up[:,n]     # Leading edge of first row of wake doublet panels
        Xw0[:,1:body['mw'][ibody2]+1]=np.matlib.repmat(np.reshape(Xw0[:,0],(m+1,1)),1,body['mw'][ibody2])+np.matlib.repmat(np.cumsum(body['wake_dx'][ibody2]),m+1,1) # Chordwise coordinates of all other wake doublet panels
        # Set y-coordinates of all spanwise wake doublet panels to corresponding
        # spanwise y-coordinates of wing panels
        Yw0=np.matlib.repmat(np.reshape(Yp_up[:,n],(m+1,1)),1,body['mw'][ibody2]+1)
        # Set z-coordinates of all wake  doublet panel vertices to corresponding
        # bound trailing edge spanwise vortex segment z-coordinates
        Zw0=np.matlib.repmat(np.reshape(Zp_up[:,n],(m+1,1)),1,body['mw'][ibody2]+1)
    else:
        Xw0=0.0
        Yw0=0.0
        Zw0=0.0
    # End if
    
    # Name of body
    name_original=body['name'][ibody2].decode('UTF-8')
    name=name_original+'_upper'  

    # Assign data to body struct array
    body['Xp0'][ibody2+1]=Xp_up
    body['Yp0'][ibody2+1]=Yp_up
    body['Zp0'][ibody2+1]=Zp_up
    body['Xc0'][ibody2+1]=Xc0
    body['Yc0'][ibody2+1]=Yc0
    body['Zc0'][ibody2+1]=Zc0 
    body['Xc0all'][ibody2+1]=Xc0all
    body['Yc0all'][ibody2+1]=Yc0all
    body['Zc0all'][ibody2+1]=Zc0all 
    body['Xw0'][ibody2+1]=Xw0
    body['Yw0'][ibody2+1]=Yw0
    body['Zw0'][ibody2+1]=Zw0
    body['mw'][ibody2+1]=body['mw'][ibody2]
    body['wake_dx'][ibody2+1]=body['wake_dx'][ibody2]
    body['nx0'][ibody2+1]=nx0
    body['ny0'][ibody2+1]=ny0
    body['nz0'][ibody2+1]=nz0
    body['tauxx0'][ibody2+1]=tauxx0
    body['tauxy0'][ibody2+1]=tauxy0
    body['tauxz0'][ibody2+1]=tauxz0
    body['tauyx0'][ibody2+1]=tauyx0
    body['tauyy0'][ibody2+1]=tauyy0
    body['tauyz0'][ibody2+1]=tauyz0
    body['nx0all'][ibody2+1]=nx0all
    body['ny0all'][ibody2+1]=ny0all
    body['nz0all'][ibody2+1]=nz0all
    body['s0'][ibody2+1]=s0
    body['s0all'][ibody2+1]=s0all
    body['m'][ibody2+1]=m
    body['n'][ibody2+1]=n
    body['c0'][ibody2+1]=body['c0'][ibody2]
    body['b'][ibody2+1]=body['b'][ibody2]
    body['yc'][ibody2+1]=body['yc'][ibody2]
    body['AR'][ibody2+1]=body['AR'][ibody2]
    body['S'][ibody2+1]=body['S'][ibody2]
    body['name'][ibody2+1]=name
    body['dir_tau'][ibody2+1]=dir_tau   
    
    return body,xyzint

def mirror_aircraft(bodyhalf,tp_body):
    # Mirrors the half-aircraft bodies in struct array bodyhalf in order to create 
    # complete (two-sided) aircraft.
    # bodyhalf: struct array containing the geometry of all the half-bodies in the flow.
    #       Each element of bodyhalf describes a different wing, fairing or 
    #       fuselage 
    # tp_body: body data type
    # body: struct array containing the geometry of all the half-bodies in the flow
    #       and their mirrors.
    #       Each element of bodyhalf describes a different wing, fairing or 
    #       fuselage 
    
    # Determine number of half-bodies
    nbody=len(bodyhalf)
    # Initialize body struct array
    body=np.zeros(2*nbody,dtype=tp_body)
    for ibody in range (0,nbody):
        # Assign right side
        body[2*ibody]=bodyhalf[ibody]
        # Create left side
        m=bodyhalf['m'][ibody]
        n=bodyhalf['n'][ibody]
        mw=bodyhalf['mw'][ibody]
        name='left'+bodyhalf['name'][ibody].decode('UTF-8')[5:]
        c0=bodyhalf['c0'][ibody]
        b=bodyhalf['b'][ibody]
        AR=bodyhalf['AR'][ibody]
        S=bodyhalf['S'][ibody]
        dir_tau=bodyhalf['dir_tau'][ibody]
        if dir == 1:    # Fuselage-type body
            Xp0=np.flipud(bodyhalf['Xp0'][ibody])
            Yp0=-np.flipud(bodyhalf['Yp0'][ibody])
            Zp0=np.flipud(bodyhalf['Zp0'][ibody])
        else:           # Wing or fin-type body
            Xp0=np.fliplr(bodyhalf['Xp0'][ibody])
            Yp0=-np.fliplr(bodyhalf['Yp0'][ibody])
            Zp0=np.fliplr(bodyhalf['Zp0'][ibody])            
        # End if
        
        # Calculate control points, normal vectors and areas
        Xc0,Yc0,Zc0,nx0,ny0,nz0,tauxx0,tauxy0,tauxz0,tauyx0,tauyy0,tauyz0,s0,cpln0,sm0,tmx0,tmy0,tmz0,sn0,tnx0,tny0,tnz0=SDPM_control_normal_tangential(Xp0,Yp0,Zp0,dir_tau)
        
        # Reshape control point coordinate matrices into vectors
        Xc0all=np.reshape(Xc0,(m*n,1),order='C')
        Yc0all=np.reshape(Yc0,(m*n,1),order='C')
        Zc0all=np.reshape(Zc0,(m*n,1),order='C')
        # Reshape normal vector component matrices into vectors
        nx0all=np.reshape(nx0,(m*n,1),order='C')
        ny0all=np.reshape(ny0,(m*n,1),order='C')
        nz0all=np.reshape(nz0,(m*n,1),order='C')
        # Reshape panel area matrix into a vector
        s0all=np.reshape(s0,(m*n,1),order='C')  
        
        # Mirror wake
        if mw != 0:
            if dir_tau == 1:    # Fuselage-type body
                Xw0=np.flipud(bodyhalf['Xw0'][ibody])
                Yw0=-np.flipud(bodyhalf['Yw0'][ibody])
                Zw0=np.flipud(bodyhalf['Zw0'][ibody])
            else: # Wing or fin-type body
                Xw0=np.fliplr(bodyhalf['Xw0'][ibody])
                Yw0=-np.fliplr(bodyhalf['Yw0'][ibody])
                Zw0=np.fliplr(bodyhalf['Zw0'][ibody])
            # End if 
        else:
            Xw0=0.0
            Yw0=0.0
            Zw0=0.0
        # End if
        # Assign data to body struct array
        body['Xp0'][2*ibody+1]=Xp0
        body['Yp0'][2*ibody+1]=Yp0
        body['Zp0'][2*ibody+1]=Zp0
        body['Xc0'][2*ibody+1]=Xc0
        body['Yc0'][2*ibody+1]=Yc0
        body['Zc0'][2*ibody+1]=Zc0 
        body['Xc0all'][2*ibody+1]=Xc0all
        body['Yc0all'][2*ibody+1]=Yc0all
        body['Zc0all'][2*ibody+1]=Zc0all 
        body['Xw0'][2*ibody+1]=Xw0
        body['Yw0'][2*ibody+1]=Yw0
        body['Zw0'][2*ibody+1]=Zw0
        body['mw'][2*ibody+1]=mw
        body['wake_dx'][2*ibody+1]=bodyhalf['wake_dx'][ibody]
        body['nx0'][2*ibody+1]=nx0
        body['ny0'][2*ibody+1]=ny0
        body['nz0'][2*ibody+1]=nz0
        body['tauxx0'][2*ibody+1]=tauxx0
        body['tauxy0'][2*ibody+1]=tauxy0
        body['tauxz0'][2*ibody+1]=tauxz0
        body['tauyx0'][2*ibody+1]=tauyx0
        body['tauyy0'][2*ibody+1]=tauyy0
        body['tauyz0'][2*ibody+1]=tauyz0
        body['nx0all'][2*ibody+1]=nx0all
        body['ny0all'][2*ibody+1]=ny0all
        body['nz0all'][2*ibody+1]=nz0all
        body['s0'][2*ibody+1]=s0
        body['s0all'][2*ibody+1]=s0all
        body['m'][2*ibody+1]=m
        body['n'][2*ibody+1]=n
        body['c0'][2*ibody+1]=c0
        body['b'][2*ibody+1]=b
        body['yc'][2*ibody+1]=-(1-bodyhalf['yc'][ibody])
        body['AR'][2*ibody+1]=AR
        body['S'][2*ibody+1]=S
        body['name'][2*ibody+1]=name
        body['dir_tau'][2*ibody+1]=dir_tau       
    # End for
    
    return body

    
    
    
 
    
    