#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Program steady_DLM_NASATM84367 calculates the steady pressure
distributions around the steady swept wing tested by Lockman and Lee
Seegmiller in NASA TM 84367, using the Doublet Lattice Method.

This code is part of the SDPMflut Python distribution.
Copyright (C) 2025 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'/Users/mbgssgd3/Documents/Python/SDPMflut_v0.72/Common/'
# 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)
from DLMgeometry import DLMgeometry_trap_fun
import DLMcalcs
# Create DLM data types
tp_trap, tp_body, _=DLMcalcs.DLMdtypes()

# Run data from wind tunnel tests of wing. Source:
# The forces and pressure distributions at subsonic speeds on a plane wing
# having 45o of sweepback, an aspect ratio of 3, and a taper ratio of 0.5,
# C. D. Kolbe and F. W. Boltz, NACA RM A51G31, 1951.
# Load experimental pressure data from .mat file
mat = scipy.io.loadmat("dataNACARMA51G31.mat")
# Total number of runs
nruns=mat['data'][0].size
# Select angle of attack values to be simulated
alpha0vec=np.array([-2,0,2,4,6,8,10,12,14,16,18])*np.pi/180.0
# Sideslip angle
beta0=0.0

# Number of bodies
nbody=1
# Initialize body struct array
body=np.zeros(nbody,dtype=tp_body)

# Input first body
ibody=0         # Index of body
name='wing'     # Name of body
# Choose numbers of panels for this wing and its wake
nhalf=20        # Number of spanwise panels per half-wing.
m=20            # Number of chordwise panels
# 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=46.67*0.0254                    # Span in m of half-wing
c0=41.47*0.0254                       # Root chord in m
lamda=0.5                       # Taper ratio
LamdaLE=48.54*np.pi/180.0;       # Sweep at leading edge in rad
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
# Calculate mean aerodynamic chord
cbarbar=2.0/3.0*c0*(lamda**2+lamda+1)/(lamda+1)
# Spanwise potiion of mean aerodynamic chord
ybarbar=bhalf/(lamda-1)*(cbarbar/c0-1.0)
# Chordwise position of leading edge of mean aerodynamic chord
xbarbar=np.tan(LamdaLE)*ybarbar
# 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='NACA64A010'
# Set airfoil parameters
teclosed=1                  # 1 for closed trailing edge, 0 otherwise
# Assemble airfoil parameter values
airfoil_params=np.array([teclosed, 0.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)

# Set characteristic chord length: root chord of wing divided by 2
bchar=trap[0]['rootchord']
# Point around which to calculate moments
xf0=xbarbar+cbarbar/4.0
yf0=0.0
zf0=0.0

# Calculate panel aspect ratio
panelAR=(c0/m)/(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=1
# Spanwise panel distribution: 1 constant, 2 denser at the wing tip(s)
linspan=1
# 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=2 # 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.
# Calculate vertices of wing panels
body=DLMgeometry_trap_fun(body,ibody,m,nhalf,mirroredwing,linchord,linspan,trap,name,rollpitchyaw,rollpitchyaw_cent,lexyz,nmin,bchar)

# 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])
# End for
axx.set_proj_type('ortho')  # FOV = 0 deg
axx.axis('equal')
axx.set_zlim(-0.02,0.02)
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 etc. for all bodies.
allbodies=DLMcalcs.allbodyindexDLM(body,bchar)

print('Calculating flutter solutions for all experimental test cases')
for irun in range (0,nruns):  
    print('Simulating run '+str(irun+1))
    
    # Set Mach number of current run
    Mach=mat['data'][0]['Mach'][irun][0][0]
    # Set Reynolds number of current run
    Reynolds=mat['data'][0]['Reynolds'][irun][0][0]
    
    # Calculate steady aerodynamic pressures and loads
    body,allbodies,barUinf,barVinf,barWinf,Cx0,Cy0,Cz0,Cl0,Cm0,Cn0=DLMcalcs.steadysolveDLM(body,allbodies,Mach,alpha0vec,beta0,xf0,yf0,zf0,install_dir)
    
    # Calculate lift coefficient and drag coefficient
    CL0=Cz0*np.cos(alpha0vec)
    CD0=Cz0*np.sin(alpha0vec)
    
    fig, axx = plt.subplots()
    axx.plot(alpha0vec*180.0/np.pi, CL0, label = "DLM")
    axx.plot(mat['data'][0]['clalpha'][irun][:,0], mat['data'][0]['clalpha'][irun][:,1], "o",label = "Exp.",)
    axx.set_xlabel("$\\alpha$")
    axx.set_ylabel("$C_L$")
    axx.grid()
    axx.legend(loc="upper left")
    plt.title("$M_{\\infty}=$" +str(Mach) + ", $Re$=" + str(Reynolds))

    fig, axx = plt.subplots()
    axx.plot(Cm0, CL0, label = "DLM")
    axx.plot(mat['data'][0]['clcm'][irun][:,0], mat['data'][0]['clcm'][irun][:,1], "o",label = "Exp.",)
    axx.set_xlabel("$C_m$")
    axx.set_ylabel("$C_L$")
    axx.grid()
    axx.legend(loc="upper left")
    plt.title("$M_{\\infty}=$" +str(Mach) + ", $Re$=" + str(Reynolds))
    
    fig, axx = plt.subplots()
    axx.plot(CD0, CL0, label = "DLM")
    axx.plot(CD0+np.min(mat['data'][0]['clcd'][irun][:,0]), CL0, label = "DLM+CD0")
    axx.plot(mat['data'][0]['clcd'][irun][:,0], mat['data'][0]['clcd'][irun][:,1], "o",label = "Exp.",)
    axx.set_xlabel("$C_D$")
    axx.set_ylabel("$C_L$")
    axx.grid()
    axx.legend(loc="upper left")
    plt.title("$M_{\\infty}=$" +str(Mach) + ", $Re$=" + str(Reynolds))
# End for
