#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Program unsteady_SDPM_LANN calculates the steady and unsteady pressure
distributions around the LANN wing in AGARD Report 702, 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
import DLMcalcs
import SDPMcalcs
# Create SDPM data types
tp_trap, tp_body, _=DLMcalcs.DLMdtypes()

## Run data from wind tunnel tests of wing. Source:
# AGARD Report 702
# Load experimental data
mat = scipy.io.loadmat("dataLANN.mat")
# Total number of runs
nruns=len(mat['data'][0])

# 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=40    # Number of spanwise panels per half-wing.
m=40        # Number of chordwise panels

# Reference chord length: root chord
cref=0.3608
# Reference span length: wing half-span
bref=1.0
# Spanwise coordinates of trapezoidal section endpoints
bhalfvec=bref*np.array([0.0, 0.2, 0.325, 0.475, 0.650, 0.825, 0.950, 1.0])
# Spanwise coordinates of measurement positions;
ymes=bref*np.array([0.2, 0.325, 0.475, 0.650, 0.825, 0.950])
# Span lengths of trapezoidal sections
bhalftrap=bhalfvec[1:]-bhalfvec[0:-1]
# Taper ratio of complete half-wing 
lambda_wing=0.4
# Root chords of trapezoidal sections
c0trap=np.interp(bhalfvec, np.array([0.0, bref]), np.array([cref, cref*lambda_wing]))
# Sweep angle at the leading edge, equal for all trapezoidal sections
LambdaLE=27.493*np.pi/180.0
# Coordinates of point through which pitch axis passes
xf0=0.224
yf0=0
zf0=0
# Airfoil section
airfoil='LANNairfoil'
# Set airfoil parameters: [Spanwise section, closed traling edge]
# Closed trailing edge
teclosed=1
# Number of trapezoidal sections
ntrap=len(bhalftrap)
# Initialize trapezoidal section struct array
trap=np.zeros(ntrap,dtype=tp_trap)

for itrap in range (0,ntrap): 
    # Define trapezoidal section itrap
    # Root chord
    c0=c0trap[itrap]
    # Half-span of trapezoidal section
    bhalf=bhalftrap[itrap]
    # Taper ratio
    lamda=c0trap[itrap+1]/c0trap[itrap]
    # Twist angle at the root
    roottwist=0.0*np.pi/180.0
    # Twist angle at the tip
    tiptwist=0.0*np.pi/180.0
    # Normalized chordwise coordinate around which wing is twisted
    twistcent=0.0
    # Dihedral angle in rad
    dihedral=0.0*np.pi/180.0  
    # Difference in leading edge root chordwise coordinate to that of the
    # previous section
    xledist=0.0    
    # Root airfoil parameters
    rootairfoilparams=np.array([itrap, teclosed])
    # Tip airfoil parameters
    tipairfoilparams=np.array([itrap+1, teclosed])
    # Arrange all data into trapezoidal sections
    trap[itrap]=np.array([(c0,xledist,bhalf,lamda,LambdaLE,roottwist,tiptwist,twistcent,dihedral,airfoil,rootairfoilparams,airfoil,tipairfoilparams)],dtype=tp_trap)  
# End for

# Calculate panel aspect ratio
panelAR=(cref/m)/(bref/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=cref/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,2):  
    print('Simulating run '+str(irun+1))
    
    # Set Mach number of current run
    Mach=mat['data'][0]['Mach'][irun][0][0]
    # Set reduced frequency
    k=mat['data'][0]['k'][irun][0][0]
    # Set mean angle of attack
    alpha0=mat['data'][0]['theta0'][irun][0]*np.pi/180.0
    
    # 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)
    # Interpolate steady pressure distribution onto spanwise measurement
    # positions
    cpspan,xspan=SDPMcalcs.cp_interpspan(body,0,body['Deltacp0'][0],ymes)

    # Plot steady pressures
    for i in range (0,len(ymes)):
        iup=np.argwhere(mat['data'][0]['Dcp'][irun][:,1] == ymes[i])
        fig, axx = plt.subplots()
        axx.plot(xspan[:,i], cpspan[:,i], label = "DLM")
        axx.plot(mat['data'][0]['Dcp'][irun][iup,0], mat['data'][0]['Dcp'][irun][iup,2], "o",label = "Exp.",)
        axx.set_xlabel("$x$")
        axx.set_ylabel("$\Delta c_p(0)$")
        axx.grid()
        axx.legend(loc="upper right")
        plt.title("$M_{\\infty}=$" +str(Mach) + ", $\\alpha_0$=" + str(alpha0*180/np.pi) + "$^o$, $2y/b$=" + f'{ymes[i]/bref:.2f}')
        plt.show()
    # End for
    
    # 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')
    # Interpolate unsteady pressure distribution onto spanwise measurement
    # positions
    cpspan,xspan=SDPMcalcs.cp_interpspan(body,0,Dcptot,ymes)

    # Plot unsteady pressures
    for i in range (0,len(ymes)):
        iup=np.argwhere(mat['data'][0]['Dcp'][irun][:,1] == ymes[i])
        fig, axx = plt.subplots(2,1)
        axx[0].plot(xspan[:,i], np.real(cpspan[:,i]), label = "DLM")
        axx[0].plot(mat['data'][0]['Dcp'][irun][iup,0], mat['data'][0]['Dcp'][irun][iup,3], "o",label = "Exp.",)
        axx[0].set_xlabel("$x$")
        axx[0].set_ylabel("$\Re(\\Delta c_p(k))/\\alpha_1$")
        axx[0].grid()
        axx[0].set_title("$M_{\\infty}=$" +str(Mach) + ", $\\alpha_0$=" + str(alpha0*180/np.pi) + "$^o$, $2y/b$=" + f'{ymes[i]/bref:.2f}')
        axx[0].legend(loc="upper right")
        axx[1].plot(xspan[:,i], np.imag(cpspan[:,i]), label = "DLM")
        axx[1].plot(mat['data'][0]['Dcp'][irun][iup,0], mat['data'][0]['Dcp'][irun][iup,4], "o",label = "Exp.",)
        axx[1].set_xlabel("$x$")
        axx[1].set_ylabel("$\Im(\\Delta c_p(k))/\\alpha_1$")
        axx[1].grid()
        plt.show()
    # End for    
# End for
