#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Program steady_SDPM_NACARML52A29 calculates the steady pressure
distributions around the wing-fuselage combination tested in 
NACA RM L52A29, using the Source and Doublet Panel Method.

This code is part of the SDPMflut Python 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/>.
"""

# Input installation directory
install_dir=r"C:\Users\pasti\Documents\Python\SDPMflut_v0.8\Common"
# install_dir=r"/Users/Username/Documents/Python/SDPMflut_v0.8/Common/" # Linux/Mac OS example
# install_dir=r"C:\Users\Username\Documents\Python\SDPMflut\Common"  # Windows example
# Import libraries and packages
import numpy as np
import matplotlib.pyplot as plt
import scipy.io
import sys
sys.path.append(install_dir)
import SDPMgeometry 
import SDPMcalcs
# Create SDPM data types
tp_trap, tp_body, _=SDPMcalcs.SDPMdtypes()

# Select wing
iwing=1     # 0: wing 45-2-0.6-006, 1: wing 45-4-0.6-006, 2: wing 45-6-0.6-006

# Run data from wind tunnel tests of wing. Source:
# Wind tunnel investigation of the aerodynamic characteristics in pitch of
# wing-fuselage combinations at high subsonic speeds, R. E. Kuhn and J. W.
# Wiggins, NACA RM L52A29
# Load data
mat = scipy.io.loadmat("dataNACARML52A9.mat")
# Total number of runs
nruns=len(mat['data'])

# Wing geometric characteristics
if iwing == 0 :      # wing 45-2-0.6-006
    Sref=3.00*0.3048**2 # Reference surface area
    bref=2.45*0.3048    # Reference span length
    cref=1.250*0.3048   # Mean aerodynamic chord
    c0=18.36*0.0254     # Root chord
    ctip=11.02*0.0254   # Tip chord
elif iwing == 1:   # wing 45-4-0.6-006
    Sref=2.25*0.3048**2 # Reference surface area
    bref=3.00*0.3048    # Reference span length
    cref=0.765*0.3048   # Mean aerodynamic chord
    c0=11.25*0.0254     # Root chord
    ctip=6.75*0.0254    # Tip chord
elif iwing == 2:   # wing 45-6-0.6-006
    Sref=2.00*0.3048**2 # Reference surface area
    bref=3.46*0.3048    # Reference span length
    cref=0.488*0.3048   # Mean aerodynamic chord
    c0=7.22*0.0254      # Root chord
    ctip=4.34*0.0254    # Tip chord
# End if

# Choose order of pressure coefficient equation
cp_order=2 # 1 for linear and 2 for second order

# Set value of mean (or steady) angle of sideslip
beta0=0.0*np.pi/180.0 # Angle of sideslip in rad

# Number of bodies
nbody=4     # Wing, wingtip, lower fuselage, upper fuselage
# Initialize bodyhalf struct array
bodyhalf=np.zeros(nbody,dtype=tp_body)

# Input first body
ibody=0         # Index of body
name='right_wing'     # Name of body
# Choose numbers of panels for this wing and its wake
nhalf=14        # Number of spanwise panels per half-wing.
mhalf=20        # Number of chordwise panels on one side of wing (upper or lower).
nchords=10      # Set length of wake in chord lengths
Dxw=4   # 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.
# Set number of trapezoidal sections for this wing
ntrap=1;
# Initialize trapezoidal section struct array
trap=np.zeros(ntrap,dtype=tp_trap)

# Input wing geometry
bhalf=bref/2.0                    # Span in m of half-wing
lamda=ctip/c0                       # Taper ratio
Lamdac4=45.0*np.pi/180.0;       # Sweep at quarter chord in rad
# Calculate the sweep angle at the leading edge
LamdaLE=np.arctan((1-lamda)*c0/4.0/bhalf+np.tan(Lamdac4))
# Mean aerodynamic chord, cmac=c0*(lambda^2+lambda+1)*2/(lambda+1)/3
cmac=c0*(lamda**2+lamda+1)*2/(lamda+1)/3
# Spanwise position of mean aerodynamic chord
ymac=(cmac-c0)*bhalf/c0/(lamda-1)
# Chordwise position of leading edge of mean aerodynamic chord
xmac=np.tan(LamdaLE)*ymac
roottwist=0.0*np.pi/180.0       # Twist angle at the root in rad
tiptwist=0.0*np.pi/180.0        # Twist angle at the root in rad
dihedral=0.0*np.pi/180.0        # Dihedral angle in rad
twistcent=0.0                   # Chordwise axis around which twist is defined
# Chordwise distance of root leading edge from previous
# trapezoidal section's tip leading edge
xledist=0.0  
# Set airfoil name (must be the filename of a function in the Common directory)
airfoil='NACA65A006'
# Assemble airfoil parameter values
airfoil_params=np.array([0.0, 0.0])
# Reference axis location
xf0=xmac+cmac/4  # Quarter chord of mean aerodynamic chord
yf0=0
zf0=0

# Arrange all data into trapezoidal sections
trap[0]=np.array([(c0,xledist,bhalf,lamda,LamdaLE,roottwist,tiptwist,twistcent,dihedral,airfoil,airfoil_params,airfoil,airfoil_params)],dtype=tp_trap)

# Calculate panel aspect ratio
panelAR=(c0/mhalf)/(bhalf/nhalf)
if panelAR < 0.1:
    sys.exit('Panel aspect ratio too low. Increase n or decrease m.')

# Minimum number of spanwise panels per trapezoidal section
nmin=3
# Chordwise panel distribution: 1 constant, 2 denser at the leading edge
linchord=0
# Spanwise panel distribution: 1 constant, 2 denser at the wing tip(s)
linspan=0
# Define root leading edge
lexyz=np.array([0, 0, 0])
# Define roll, pitch and yaw angles
rollpitchyaw=np.array([0, 0, 0])*np.pi/180;
# Define roll, pitch and yaw centre (x,y,z position of rotation centre)
rollpitchyaw_cent=np.array([0, 0, 0]);
# Input body description
mirroredwing=1 # 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.
# dir_tau is the direction in which the unit tangent vector  for this wing
# (tauxx, tauxy, tauxz) has a zero component
dir_tau=2
# Calculate vertices of wing panels
bodyhalf=SDPMgeometry.SDPMgeometry_trap_fun(bodyhalf,ibody,mhalf,nhalf,nchords,Dxw,mirroredwing,linchord,linspan,trap,name,dir_tau,rollpitchyaw,rollpitchyaw_cent,lexyz,nmin,cref)
# Calculate vertices of wingtip panels
bodyhalf=SDPMgeometry.makewingtips(bodyhalf,ibody,ibody,mirroredwing)

# Create right fuselage
# Distance of fuselage nose quarter of wing mean aerodynamic chord
nose_x=30*0.0254
# Load geometry file
matfus = scipy.io.loadmat("NACARML5A29_fuselage_sections.mat")
fuselage_sections=matfus['fuselage_sections'][0]
num_sec=len(fuselage_sections)-1
# Lengthwise positions of trapezoidal sections
xfus=np.zeros(num_sec+1)
for i in range (0,num_sec+1):
    xfus[i]=fuselage_sections['x'][i][0][0]
# End for

# Set up fuselage
ibody=2
# Name of body
name='right_fuselage';
# Set numbers of lengthwise and heightwise panels for fuselage
m=20 # Number of lengthwise panels (nose to tail)
n=20 # Number of heightwise panels (bottom to top)
nlengths=5 # Length of wake in fuselage lengths
Dxw=0.0 # No wake
# Minimum number of spanwise panels per trapezoidal section
nmin=1 
# Lengthwise panel distribution: 1 constant, 2 denser at the leading edge
linlength=0
# Distance of nose leading edge from (0,0,0)
nosexyz=np.array((-nose_x+xmac+cmac/4, 0.0, 0.0))
# Calculate vertices of fuselage panels
bodyhalf=SDPMgeometry.SDPMgeometry_fuselage(bodyhalf,ibody,m,nlengths,Dxw,n,nmin,linlength,xfus,fuselage_sections,nosexyz,name)

# Determine intersections between wing and fuselage
ibody1=0 # Index of wing in body
ibody2=2 # Index of fuselage in body
j1min=0  # Spanwise section of the wing at which to start looking for intersections.
bodyhalf,xyzint=SDPMgeometry.wing_fuselage_intersections(bodyhalf,ibody1,ibody2,j1min,install_dir)

# Mirror all body halves to obtain a complete aircraft
body=SDPMgeometry.mirror_aircraft(bodyhalf,tp_body)

# Plot all bodies
fig, axx = plt.subplots(subplot_kw={"projection": "3d"})
for i in range (0,len(body)):
    axx.plot_surface(body['Xp0'][i], body['Yp0'][i], body['Zp0'][i],edgecolor='royalblue', lw=0.5,alpha=0.2)
    if body['mw'][i] > 0:
        axx.plot_surface(body['Xw0'][i], body['Yw0'][i], body['Zw0'][i],edgecolor='red', lw=0.5,alpha=0.2, rstride=1, cstride=1)
    # End if    
# End for
axx.set_proj_type('ortho')  # FOV = 0 deg
axx.axis('equal')
axx.set_xlabel("$x$", labelpad=10)
axx.set_ylabel("$y$", labelpad=10)
axx.set_zlabel("$z$", labelpad=-1)
axx.view_init(26, -120)
plt.show()

# Assemble the indices of the body panels, spanwise body panels, wake
# panels etc. for all bodies.
allbodies=SDPMcalcs.allbodyindex(body)

alpha0=np.arange(-6,13,2)*np.pi/180.0
dCmdCL=np.zeros(nruns)

for irun in range (0,nruns):
    print('Simulating run '+str(irun+1))
    
    # Set Mach number of current run
    Mach=mat['data']['Mach'][irun][iwing][0][0]
    # Calculate subsonic compressibility factor
    beta=np.sqrt(1-Mach**2) 
    
    # Calculate steady aerodynamic pressures and loads
    body,allbodies,Aphi,Bphi,Cphi,barUinf,barVinf,barWinf,Cx0,Cy0,Cz0,Cl0,Cm0,Cn0=SDPMcalcs.steadysolve(body,allbodies,cp_order,Mach,beta,alpha0,beta0,xf0,yf0,zf0,Sref,cref,bref,install_dir)
    # Calculate lift and drag
    CL0=-Cx0*np.sin(alpha0)+Cz0*np.cos(alpha0)
    CD0=Cx0*np.cos(alpha0)+Cz0*np.sin(alpha0)
    
    fig, axx = plt.subplots()
    axx.plot(alpha0*180.0/np.pi,CL0, label = "SDPM")
    axx.plot(mat['data']['CL'][irun][iwing][:,0],mat['data']['CL'][irun][iwing][:,1], "o", label = "Exp.")
    axx.set_xlabel("$\\alpha$")
    axx.set_ylabel("$C_L$")
    axx.grid()
    axx.legend(loc="lower right")
    plt.title("Wing No. " +str(iwing) +", $M_{\\infty}=$" +str(Mach))
    plt.show()    
    
    fig, axx = plt.subplots()
    axx.plot(CL0,Cm0, label = "SDPM")
    axx.plot(mat['data']['CM'][irun][iwing][:,0],mat['data']['CM'][irun][iwing][:,1], "o", label = "Exp.")
    axx.set_xlabel("$C_L$")
    axx.set_ylabel("$C_M$")
    axx.grid()
    axx.legend(loc="lower right")
    plt.title("Wing No. " +str(iwing) +", $M_{\\infty}=$" +str(Mach))
    plt.show()   
    
    fig, axx = plt.subplots()
    axx.plot(CL0,CD0, label = "SDPM $C_D+C_{D_0}$")
    axx.plot(CL0,CD0+mat['data']['CD'][irun][iwing][1,1], label = "SDPM")
    axx.plot(mat['data']['CD'][irun][iwing][:,0],mat['data']['CD'][irun][iwing][:,1], "o", label = "Exp.")
    axx.set_xlabel("$C_L$")
    axx.set_ylabel("$C_D$")
    axx.grid()
    axx.legend(loc="upper left")
    plt.title("Wing No. " +str(iwing) +", $M_{\\infty}=$" +str(Mach))
    plt.show()      
# End for

    