#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Program unsteady_DLM_NLR7301 calculates the steady and unsteady pressure
distributions around the 2D wing undergoing pitching or control surface 
oscillations, as described in report AGARD-R-702 (Data set 4), 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"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)
from DLMgeometry import DLMgeometry_trap_fun
from SDPMcalcs import modeshape_assemble
import DLMcalcs
# Create SDPM data types
tp_trap, tp_body, _=DLMcalcs.DLMdtypes()

# Run data from wind tunnel tests of wing. Source:
# Zwaan RJ 1982 Data set 4: NLR 7301 supercritical airfoil. oscillatory 
# pitching and oscillating flap. In Compendium of unsteady aerodynamic 
# measurements (ed. Various) number AGARD-R-702 AGARD pp. 4–1 – 4–25. Only 
# pitching cases are included, an oscillating control surface case will be
# included at a later date.
# Free stream Mach number
Machdata=np.array([0.499, 0.499, 0.498, 0.696, 0.696]) # Cases 4 and 5 are supercritical and feature shock waves on the surface
# Mean angle of attack in rad
alpha0data=np.array([0.85, 0.85, 0.85, 3.00, 3.00])*np.pi/180.0
# Pitch amplitude in rad
alpha1data=np.array([0.50, 0.55, 0.44, 0.50, 0.42])*np.pi/180.0
# Reduced frequency
kdata=np.array([0.000, 0.098, 0.262, 0.000, 0.072])
# Total number of runs
nruns=Machdata.size
# Load experimental pressure data from .mat file
mat = scipy.io.loadmat("dataNLR7301.mat")

# 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=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=50    # 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=18*0.0254*10              # Span in m of half-wing
c0=18*0.0254                    # Root chord in m
lamda=1.0                       # Taper ratio
LamdaLE=0.0*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
# 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='NLR7301'
# Set airfoil parameters
dataset=2   # Select the measured airfoil data given in AGARD-R-702.
teclosed=1  # Closed trailing edge
# Assemble airfoil parameter values
airfoil_params=np.array([dataset, teclosed])
# Pitch axis location
xf0=0.4*c0 
yf0=0.0; 
zf0=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)

# 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.')

# 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]);
# 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])
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)

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=Machdata[irun]
    # Set reduced frequency
    k=kdata[irun]
    # Set mean angle of attack
    alpha0=np.array([alpha0data[irun]])
    # Set pitch amplitude
    alpha1=np.array([alpha1data[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,0.0,0.0,0.0,install_dir)

    DeltaCP_exp=mat['data']['CP_lo'][0][irun]-mat['data']['CP_up'][0][irun]

    # Plot steady pressure
    fig, axx = plt.subplots()
    axx.plot(body['Xc0'][0][:,body['n'][0]//2]/c0, body['Deltacp0'][0][:,body['n'][0]//2], label = "DLM")    # Plot experimental pressure measurements
    axx.plot(np.ndarray.flatten(mat['data']['x'][0][irun]), DeltaCP_exp, "o", label = "Exp. upper")    # Plot experimental pressure measurements
    axx.set_xlabel("$x/c_0$", labelpad=10)
    axx.set_ylabel("$\Delta c_{p_0}$", labelpad=10)
    axx.grid()
    axx.legend(loc="upper right")
    plt.title("$M_{\\infty}=$" +str(Mach) + ", $\\alpha_0$=" + str(alpha0[0]*180/np.pi) + "$^o$")
    plt.show()
    
    # Calculate the unsteady pressure coefficients
    Deltacphdot,Deltacpalpha,Deltacpalphadot=DLMcalcs.unsteadysolve_pitchplungeDLM(body,allbodies,k,Mach,xf0,yf0,install_dir)
    # Calculate total oscillatory pressure jump distribution
    Dcptot=(Deltacpalpha+1j*k*Deltacpalphadot)
    Dcptot=np.reshape(Dcptot,(body['m'][0],body['n'][0]),order='C')

    DeltaCPRE_exp=mat['data']['CPRE_lo'][0][irun]-mat['data']['CPRE_up'][0][irun]
    DeltaCPIM_exp=mat['data']['CPIM_lo'][0][irun]-mat['data']['CPIM_up'][0][irun]  

    # Plot real part of unsteady pressure
    fig, axx = plt.subplots()
    axx.plot(body['Xc0'][0][:,body['n'][0]//2]/c0, np.real(Dcptot[:,body['n'][0]//2]), label = "DLM")    # Plot experimental pressure measurements
    axx.plot(np.ndarray.flatten(mat['data']['x'][0][irun]), DeltaCPRE_exp, "o", label = "Exp. upper")    # Plot experimental pressure measurements
    axx.set_xlabel("$x/c_0$", labelpad=10)
    axx.set_ylabel("$\Re(\Delta c_{p}(k))$", labelpad=10)
    axx.grid()
    axx.legend(loc="upper right")
    plt.title("$M_{\\infty}=$" +str(Mach) + ", $\\alpha_0$=" + str(alpha0[0]*180/np.pi) + "$^o$, $k=$" + str(k))
    plt.show()

    # Plot imaginary part of unsteady pressure
    fig, axx = plt.subplots()
    axx.plot(body['Xc0'][0][:,body['n'][0]//2]/c0, np.imag(Dcptot[:,body['n'][0]//2]), label = "DLM")    # Plot experimental pressure measurements
    axx.plot(np.ndarray.flatten(mat['data']['x'][0][irun]), DeltaCPIM_exp, "o", label = "Exp. upper")    # Plot experimental pressure measurements
    axx.set_xlabel("$x/c_0$", labelpad=10)
    axx.set_ylabel("$\Im(\Delta c_{p}(k))$", labelpad=10)
    axx.grid()
    axx.legend(loc="lower right")
    plt.title("$M_{\\infty}=$" +str(Mach) + ", $\\alpha_0$=" + str(alpha0[0]*180/np.pi) + "$^o$, $k=$" + str(k))
    plt.show()     
# End for
