#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Program flutter_DLM_NACA0012 calculates the flutter boundary of the NASA
Benchmark Supercritical wing with pitch and plunge degrees of freedom, using 
the Doublet Lattice 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 sys
sys.path.append(install_dir)
from DLMgeometry import DLMgeometry_trap_fun
import flutsol
import DLMcalcs
# Create DLM data types
tp_trap, tp_body, _=DLMcalcs.DLMdtypes()

# Choose tests in air or R-12
R12=1  # 0 for air, 1 for R-12
                    
# Run data from wind tunnel flutter tests. Source:
# Experimental unsteady pressures at flutter on the supercritical wing
# benchmark model. Bryan E. Dansbeny et al, AIAA-93-1592-CP, pp. 2504-2514.
# Test Cases for Flutter of the Benchmark Models Rectangular Wings on the Pitch
# and Plunge Apparatus. Robert M. Bennett. Defense Technical Information Center
# Compilation Part Notice ADPO10713.
# Chordwise reference length in m
cref=16*0.0254 
if R12 == 0:
    # Free stream Mach number
    Machdata=np.array([0.319, 0.509, 0.730, 0.769, 0.326, 0.513, 0.725, 0.766, 0.335, 0.503, 0.619, 0.679, 0.738, 0.762])
    # Mean pitch angle
    alpha0data=np.array([-0.1, 0.0, -0.2, 0.0, 1.0, 1.0, 1.2, 1.2, 0.0, -0.2, -0.4, -0.1, -0.1, -0.1])*np.pi/180.0
    # Speed of sound data
    asounddata=np.array([1142.2, 1129.6, 1120.2, 1115.6, 1140.2, 1130.4, 1121.1, 1117.8, 1139, 1132, 1125, 1123, 1118, 1116])*0.3048
    # Free stream air density (kg/m^3)
    rhodata=515.37882*np.array([2.113, 0.922, 0.516, 0.459, 2.097, 0.930, 0.516, 0.461, 2.132, 0.987, 0.693, 0.588, 0.522, 0.480])*1e-3
    # Mass ratio Mh./(pi*rhodata*b*c0^2/4)
    mubardata=np.array([776, 1779, 3178, 3569, 782, 1763, 3180, 3556, 769, 1660, 2366, 2789, 3142, 3413])*1.0
    # Flutter speed (m/s)
    Uflutdata=np.array([364.4, 574.7, 817.6, 857.6, 372.0, 580.0, 812.4, 856.8, 381.9, 569.5, 696.9, 763.0, 825.0, 850.8])*0.3048
    # Flutter frequency (rad/s)
    freqflutdata=np.array([4.60, 4.47, 4.23, 4.14, 4.55, 4.43, 4.19, 4.09, 4.55, 4.47, 4.39, 4.36, 4.27, 4.23])*2.0*np.pi
    # Frequency of first torsion mode from experiment
    omega_alpha=5.25*2.0*np.pi 
    # Nondimensional flutter speed Uflutdata./(c0/2*34.2*sqrt(mubardata)). 
    # Different values in different sources and they both do not agree with the definition of Ustar
    # Ustardata=[0.595 0.620 0.622 0.692 0.605 0.628 0.655 0.653 0.626 0.636 0.652 0.657 0.669 0.662]; According to Rivera
    # Ustardata=[0.68 0.73 0.83 0.82 0.70 0.74 0.82 0.84 0.72 0.75 0.78 0.79 0.82 0.82]; According to AIAA-93-1592-CP
    # Recalculate experimental Ustar values
    Ustardata=np.divide(Uflutdata,cref/2.0*omega_alpha*np.sqrt(mubardata))
    # Flutter frequency ratio freqflutdata/wn(2)
    freqratdata=np.array([0.876, 0.851, 0.806, 0.789, 0.867, 0.844, 0.798, 0.779,  0.867, 0.851, 0.836, 0.830, 0.813, 0.806])
    # Flutter dynamic pressure data (N/m^2)
    dynpressdata=np.array([140.3, 152.2, 172.4, 168.9, 145.1, 156.4, 170.2, 169.3, 155.5, 160.1, 168.2, 171.1, 177.6, 173.8])*4.44822/0.3048**2
    # Static pressure calculated from speed of sound and density
    pressdata=asounddata**2.*rhodata/1.4
    # Reduced frequency
    kflutdata=np.array([0.0529, 0.0326, 0.0217, 0.0202, 0.0513, 0.0320, 0.0216, 0.0200, 0.0499, 0.0329, 0.0264, 0.0239, 0.0217, 0.0208])
else:
    # Free stream Mach number
    Machdata=np.array([0.435, 0.579, 0.689, 0.742, 0.797, 0.803, 0.799, 0.817, 0.823])
    # Mean pitch angle
    alpha0data=np.array([-0.3, -0.2, -0.2, -0.1, 0.0, 0.0, 1.3, 0.0, 0.0])*np.pi/180.0
    # Speed of sound data
    asounddata=np.array([506.5, 508.0, 507.9, 506.6, 510.9, 510.0, 509.7, 511.5, 510.0])*0.3048
    # Free stream air density (kg/m^3)
    rhodata=515.37882*np.array([6.482, 3.752, 2.764, 2.392, 2.080, 2.038, 2.012, 2.050, 2.016])*1e-3
    # Mass ratio Mh./(pi*rhodata*b*c0^2/4)
    mubardata=np.array([253, 437, 593, 685, 788, 805, 815, 800, 813])*1.0
    # Flutter speed (m/s)
    Uflutdata=np.array([220.4, 294.2, 350.0, 375.7, 407.4, 409.3, 407.3, 417.7, 419.7])*0.3048
    # Nondimensional flutter speed Uflutdata./(c0/2*34.2*sqrt(mubardata))
    Ustardata=np.array([0.630, 0.640, 0.654, 0.653, 0.660, 0.656, 0.649, 0.672, 0.669])
    # Flutter frequency (rad/s)
    freqflutdata=np.array([4.53, 4.45, 4.35, 4.30, 4.14, 4.15, 4.14, 4.05, 4.13])*2.0*np.pi
    # Frequency of first torsion mode from experiment
    omega_alpha=5.25*2.0*np.pi 
    # Flutter frequency ratio freqflutdata/wn(2)
    freqratdata=np.array([0.863, 0.848, 0.829, 0.819, 0.789, 0.790, 0.789, 0.771, 0.787])
    # Flutter dynamic pressure data (N/m^2)
    dynpressdata=np.array([157.4, 162.4, 169.3, 168.8, 172.6, 170.7, 166.9, 178.8, 177.5])*4.44822/0.3048**2
    # Static pressure calculated from speed of sound and density
    pressdata=asounddata**2.*rhodata/1.4
    # Reduced frequency
    kflutdata=np.array([0.0861, 0.0633, 0.0521, 0.0480, 0.0425, 0.0424, 0.0425, 0.0407, 0.0412])
# End if
# Sideslip angle
beta0=0.0
# Total number of runs
nruns=len(Machdata);

# Structural dynamic properties measured from experiment
mass=6.1*14.5939 # Mass of wing in kg
Ialpha=2.7*14.5939*0.3048**2 # Moment of inertia around pitch axis in kg*m^2
fn=np.array([3.32, 5.25])   # Natural frequencies in plunge and pitch in Hz
wn=2*np.pi*fn             #  Natural frequencies in plunge and pitch in rad/s
zeta0=np.array([0.000, 0.000]) # Damping ratios in plunge and pitch
Kh=2637*4.44822/0.3048; # Plunge stiffness in N/m
Kalpha=2964*0.3048*4.44822; # Pitch stiffness in Nm/rad
# Number of structural modes
nmodes=2

# Select reduced frequency values
kvec=np.array([0.0, 0.001, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.1])

# Select airspeed range in m/s
Uv=np.linspace(50,320,num=101)

# 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
# Number of trapezoidal sections
ntrap=1;
# Initialize trapezoidal section struct array
trap=np.zeros(ntrap,dtype=tp_trap)

# Input wing geometry
bhalf=32*0.0254                 # Span in m of half-wing
c0=16*0.0254                    # Root chord in m
lamda=1.0                       # Taper ratio
LamdaLE=0.0*np.pi/180.0;        # Sweep at leading edge
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  
# Coordinates of points through which pitch axis passes
xf0=c0/2.0
yf0=0.0
zf0=0.0
# Set airfoil name (must be the filename of a function in the Common directory)
airfoil='NASASC20414'
# Set airfoil parameters
teclosed=1                  # 1 for closed trailing edge, 0 otherwise
# Assemble airfoil parameter values
airfoil_params=np.array([teclosed, 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.')

# 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']/2.0
# 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])
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)

# The pitch-plunge wing has two degrees of freedom and therefore two modes of 
# vibration
nmodes=2
# Parameter to determine if the structural model concerns a half wing or a
# full wing.
halfwing=1 # halfwing=1: half-wing. halfwing=0: full wing

# 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_xlabel("$x$", labelpad=10)
axx.set_ylabel("$y$", labelpad=10)
axx.set_zlabel("$z$", labelpad=-2)
axx.view_init(26, -120)
plt.show()

# Assemble the indices of the body panels, spanwise body panels, wake
# panels etc. for all bodies.
allbodies=DLMcalcs.allbodyindexDLM(body,bchar)

# Calculate structural matrices
A=np.matrix([[mass, 0.0],[0.0, Ialpha]]) # Structural mass matrix
E=np.matrix([[Kh, 0.0],[0.0, Kalpha]])   # Structural stiffness matrix
C=np.matrix([[2*mass*wn[0]*zeta0[0], 0.0],[0.0, 2*Ialpha*wn[1]*zeta0[1]]])   # Structural damping matrix

# Initialize results arrays for all runs
Ustarvec=np.zeros((nruns))    # Flutter speed index
Uflutvec=np.zeros((nruns))    # Flutter speed in m/s
dynpressvec=np.zeros((nruns)) # Flutter dynamic pressure in Pa
freqflutvec=np.zeros((nruns)) # Flutter frequency in rad/s
kflutvec=np.zeros((nruns))    # Reduced flutter frequency

print('Calculating flutter solutions for all experimental test cases')
for irun in range (0,nruns):
    print('')
    print('Simulating run '+str(irun+1))

    # Set Mach number of current run
    Mach=Machdata[irun]
    # Set mass ratio of current run
    mubar=mubardata[irun]
    # Set density of current run
    rho=rhodata[irun]
    # Set mean angle of attack
    alpha0=np.array([alpha0data[irun]])
    
    # Calculate steady aerodynamic pressures and loads
    body,allbodies,barUinf,barVinf,barWinf,Cx0,Cy0,Cz0,Cl0,Cm0,Cn0=DLMcalcs.steadysolveDLM(body,allbodies,Mach,alpha0,beta0,xf0,yf0,zf0,install_dir)
        
    # Calculate flutter solution for pitch-plunge motion
    Uflut,freqflut,kflut,dynpressflut,omega,zeta=flutsol.flutsolve_pitchplungeDLM(body,allbodies,kvec,Uv,nmodes,Mach,A,C,E,rho,wn,halfwing,xf0,yf0,zf0,install_dir)
    
    # Print out flutter solution
    if Uflut != 0: # If there is a flutter point
        # Store flutter data for this run
        Uflutvec[irun]=Uflut
        Ustarvec[irun]=Uflut/(c0/2.0*omega_alpha*np.sqrt(mubar))
        dynpressvec[irun]=dynpressflut
        freqflutvec[irun]=freqflut
        kflutvec[irun]=kflut
        # Compare to experimental data        
        print('Flutter speed (m/s)')
        print('DLM      Exp.')
        print(np.array([Uflut,Uflutdata[irun]]))
        print('Flutter frequency (rad/s)')
        print('DLM      Exp.')
        print(np.array([freqflut,freqflutdata[irun]]))
    else:
        print('Could not find a flutter point for run '+str(irun))
    # End if 
# End loop for irun

if R12 == 0:
    # Plot flutter airspeed
    fig, axx = plt.subplots()
    axx.plot(Machdata[0:4], Uflutvec[0:4], label = "DLM $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata[4:8], Uflutvec[4:8], label = "DLM $\\alpha_0=1^{\\circ}$")
    axx.plot(Machdata[8:nruns], Uflutvec[8:nruns], label = "DLM grit")
    axx.plot(Machdata[0:4], Uflutdata[0:4], "o", label = "Exp. $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata[4:8], Uflutdata[4:8], "o", label = "Exp. $\\alpha_0=1^{\\circ}$")
    axx.plot(Machdata[8:nruns], Uflutdata[8:nruns], "o", label = "Exp. grit")
    axx.set_xlabel("$M_{\\infty}$")
    axx.set_ylabel("$U_F$ (m/s)")
    axx.grid()
    axx.legend(loc="upper left")
    
    # Plot flutter frequency
    fig, axx = plt.subplots()
    axx.plot(Machdata[0:4], freqflutvec[0:4], label = "DLM $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata[4:8], freqflutvec[4:8], label = "DLM $\\alpha_0=1^{\\circ}$")
    axx.plot(Machdata[8:nruns], freqflutvec[8:nruns], label = "SDPM grit")
    axx.plot(Machdata[0:4], freqflutdata[0:4], "o", label = "Exp. $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata[4:8], freqflutdata[4:8], "o", label = "Exp. $\\alpha_0=1^{\\circ}$")
    axx.plot(Machdata[8:nruns], freqflutdata[8:nruns], "o", label = "Exp. grit")
    axx.set_xlabel("$M_{\\infty}$")
    axx.set_ylabel("$\\omega_F$ (rad/s)")
    axx.grid()
    axx.legend(loc="lower left")
    
    # Plot flutter speed index
    fig, axx = plt.subplots()
    axx.plot(Machdata[0:4], Ustarvec[0:4], label = "DLM $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata[4:8], Ustarvec[4:8], label = "DLM $\\alpha_0=1^{\\circ}$")
    axx.plot(Machdata[8:nruns], Ustarvec[8:nruns], label = "DLM grit")
    axx.plot(Machdata[0:4], Ustardata[0:4], "o", label = "Exp. $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata[4:8], Ustardata[4:8], "o", label = "Exp. $\\alpha_0=1^{\\circ}$")
    axx.plot(Machdata[8:nruns], Ustardata[8:nruns], "o", label = "Exp. grit")
    axx.set_xlabel("$M_{\\infty}$")
    axx.set_ylabel("$Q^*_F$")
    axx.grid()
    axx.legend(loc="upper left")
    
    # Plot flutter frequency ratio
    fig, axx = plt.subplots()
    axx.plot(Machdata[0:4], freqflutvec[0:4]/omega_alpha, label = "DLM $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata[4:8], freqflutvec[4:8]/omega_alpha, label = "DLM $\\alpha_0=1^{\\circ}$")
    axx.plot(Machdata[8:nruns], freqflutvec[8:nruns]/omega_alpha, label = "DLM grit")
    axx.plot(Machdata[0:4], freqratdata[0:4], "o", label = "Exp. $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata[4:8], freqratdata[4:8], "o", label = "Exp. $\\alpha_0=1^{\\circ}$")
    axx.plot(Machdata[8:nruns], freqratdata[8:nruns], "o", label = "Exp. grit")
    axx.set_xlabel("$M_{\\infty}$")
    axx.set_ylabel("$\\omega_F/\\omega_alpha$")
    axx.grid()
    axx.legend(loc="lower left")
    
    # Plot flutter dynamic pressure
    fig, axx = plt.subplots()
    axx.plot(Machdata[0:4], dynpressvec[0:4], label = "DLM $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata[4:8], dynpressvec[4:8], label = "DLM $\\alpha_0=1^{\\circ}$")
    axx.plot(Machdata[8:nruns], dynpressvec[8:nruns], label = "DLM grit")
    axx.plot(Machdata[0:4], dynpressdata[0:4], "o", label = "Exp. $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata[4:8], dynpressdata[4:8], "o", label = "Exp. $\\alpha_0=1^{\\circ}$")
    axx.plot(Machdata[8:nruns], dynpressdata[8:nruns], "o", label = "Exp. grit")
    axx.set_xlabel("$M_{\\infty}$")
    axx.set_ylabel("$1/2\\rho Q_F^2$ (Pa)")
    axx.grid()
    axx.legend(loc="upper left")
    
    # Plot flutter reduced frequency
    fig, axx = plt.subplots()
    axx.plot(Machdata[0:4], kflutvec[0:4], label = "DLM $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata[4:8], kflutvec[4:8], label = "DLM $\\alpha_0=1^{\\circ}$")
    axx.plot(Machdata[8:nruns], kflutvec[8:nruns], label = "DLM grit")
    axx.plot(Machdata[0:4], kflutdata[0:4], "o", label = "Exp. $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata[4:8], kflutdata[4:8], "o", label = "Exp. $\\alpha_0=1^{\\circ}$")
    axx.plot(Machdata[8:nruns], kflutdata[8:nruns], "o", label = "Exp. grit")
    axx.set_xlabel("$M_{\\infty}$")
    axx.set_ylabel("$k$")
    axx.grid()
    axx.legend(loc="upper right")
else:
    # Plot flutter airspeed
    fig, axx = plt.subplots()
    axx.plot(Machdata, Uflutvec, label = "DLM $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata, Uflutdata, "o", label = "Exp. $\\alpha_0=0^{\\circ}$")
    axx.set_xlabel("$M_{\\infty}$")
    axx.set_ylabel("$U_F$ (m/s)")
    axx.grid()
    axx.legend(loc="upper left")

    # Plot flutter frequency
    fig, axx = plt.subplots()
    axx.plot(Machdata, freqflutvec, label = "DLM $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata, freqflutdata, "o", label = "Exp. $\\alpha_0=0^{\\circ}$")
    axx.set_xlabel("$M_{\\infty}$")
    axx.set_ylabel("$\\omega_F$ (rad/s)")
    axx.grid()
    axx.legend(loc="lower left")

    # Plot flutter speed index
    fig, axx = plt.subplots()
    axx.plot(Machdata, Ustarvec, label = "DLM $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata, Ustardata, "o", label = "Exp. $\\alpha_0=0^{\\circ}$")
    axx.set_xlabel("$M_{\\infty}$")
    axx.set_ylabel("$Q^*_F$")
    axx.grid()
    axx.legend(loc="upper left")

    # Plot flutter frequency ratio
    fig, axx = plt.subplots()
    axx.plot(Machdata, freqflutvec/wn[1], label = "DLM $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata, freqratdata, "o", label = "Exp. $\\alpha_0=0^{\\circ}$")
    axx.set_xlabel("$M_{\\infty}$")
    axx.set_ylabel("$\\omega_F/\\omega_alpha$")
    axx.grid()
    axx.legend(loc="lower left")

    # Plot flutter dynamic pressure
    fig, axx = plt.subplots()
    axx.plot(Machdata, dynpressvec, label = "DLM $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata, dynpressdata, "o", label = "Exp. $\\alpha_0=0^{\\circ}$")
    axx.set_xlabel("$M_{\\infty}$")
    axx.set_ylabel("$1/2\\rho Q_F^2$ (Pa)")
    axx.grid()
    axx.legend(loc="upper left")

    # Plot flutter reduced frequency
    fig, axx = plt.subplots()
    axx.plot(Machdata, kflutvec, label = "DLM $\\alpha_0=0^{\\circ}$")
    axx.plot(Machdata, kflutdata, "o", label = "Exp. $\\alpha_0=0^{\\circ}$")
    axx.set_xlabel("$M_{\\infty}$")
    axx.set_ylabel("$k$")
    axx.grid()
    axx.legend(loc="upper right")
# End if
