#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Program unsteady_SDPM_Clipped_Delta calculates the steady and unsteady 
pressure distributions around the clipped Delta wing wing undergoing 
pitch oscillations tested by Hess et al in NASA-TP-2594, using the Source 
and Doublet Panel Method.

This code is part of the SDPMflut Python distribution.
Copyright (C) 2026 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 SDPMgeometry import SDPMgeometry_trap_fun
from SDPMgeometry import makewingtips
import SDPMcalcs
# Create SDPM data types
tp_trap, tp_body, _=SDPMcalcs.SDPMdtypes()

# Run data from wind tunnel tests of wing. Source:
# R. W. Hess and F. W. Cazier, Jr. and E. C. Wynne, Steady and Unsteady
# Transonic Pressure Measurements on a Clipped Delta Wing for Pitching and
# Control Surface Oscillations, NASA-TP-2594, 1988
# and
# R. M. Bennett and C. E. Walker, Computational Test Cases for a Clipped
# Delta Wing With Pitching and Trailing-Edge Control Surface Oscillations, 
# NASA/TM-1999-209104, 1999
# Free stream Mach number
Machdata=np.array([0.4030]) # The other test cases are highly transonic
# Mean angle of attack in rad
alpha0data=np.array([3.04])*np.pi/180
# Pitch amplitude in rad
alpha1data=np.array([0.46])*np.pi/180
# Frequency in rad/s
freqdata=np.array([7.99])*2*np.pi
# Reduced frequency based on mean chord
kdata=np.array([0.3830])
# Dynamic pressure in N
dynpressdata=np.array([47.8803])*78.3
# Reynolds number
Reynoldsdata=np.array([9.22e6])
# Steady spanwise measurement positions
ysmes=np.array([0.380746, 0.619506, 0.67183, 0.794258, 0.974852])
# Unsteady spanwise measurement positions
yumes=np.array([0.385826, 0.624586, 0.67564, 0.799592, 0.979678])

# Total number of runs
nruns=Machdata.size
# Load experimental pressure data from .mat file
mat = scipy.io.loadmat("data_Clipped_Delta.mat")

# 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=3     # One for the wing and two for its wingtips
# 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=10        # 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=1   # 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.
# Number of trapezoidal sections
ntrap=1;
# Initialize trapezoidal section struct array
trap=np.zeros(ntrap,dtype=tp_trap)

# Input wing geometry
bhalf=45.08*0.0254              # Span in m of half-wing
c0=63.55*0.0254                    # Root chord in m
lamda=9.03/63.55                      # Taper ratio
LamdaLE=50.4*np.pi/180.0;        # Sweep at leading edge in rad
# Geometric mean chord
cbar=(c0+c0*lamda)/2.0 
# Mean aerodynamic chord, cmac=c0*(lambda^2+lambda+1)*2/(lambda+1)/3
cmac=c0*(lamda**2.0+lamda+1.0)*2.0/(lamda+1.0)/3.0
# Spanwise position of mean aerodynamic chord
ymac=(cmac-c0)*bhalf/c0/(lamda-1.0)
# 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='biconvex'
# Set airfoil parameters
thick=0.06  # Airfoil thickness to chord ratio
# Assemble airfoil parameter values
airfoil_params=np.array([thick, 0.0])
# Reference areas and lengths for all bodies in the flow
Sref=c0*(1.0+lamda)*bhalf   # Reference surface area
bref=2*bhalf;               # Reference span length
cref=cbar;                  # This test case uses the geometric mean chord, see NASA-TP-2594
# Pitch axis location
xf0=65.2/100*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)

# Mean chord
cbar=(c0+c0*lamda)/2

# 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=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.
# 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
body=SDPMgeometry_trap_fun(body,ibody,mhalf,nhalf,nchords,Dxw,mirroredwing,linchord,linspan,trap,name,dir_tau,rollpitchyaw,rollpitchyaw_cent,lexyz,nmin,cref)
# Calculate vertices of wingtip panels
body=makewingtips(body,ibody,ibody,mirroredwing)

# 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=SDPMcalcs.allbodyindex(body)

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]
    # Calculate reduced frequency based on the root chord. The way the SDPM
    # is coded, the reduced frequency must always be based on the root chord.
    k=kdata[irun]*c0/cbar
    # Set mean angle of attack
    alpha0=np.array([alpha0data[irun]])
    # Set pitch amplitude
    alpha1=np.array([alpha1data[irun]])
    # 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)

    # Interpolate steady pressure distribution onto spanwise measurement
    # positions
    [cpspan,xspan]=SDPMcalcs.cp_interpspan(body,0,body['cp0'][0],ysmes)

    # Plot steady pressures
    for i in range (0,len(ysmes)):
        iup=np.argwhere(np.round(mat['ysall']*10000.0) == np.round(ysmes[i]*10000.0))
        iup=iup[:,0]
        fig, axx = plt.subplots()
        axx.plot(xspan[:,i], cpspan[:,i], label = "SDPM")
        axx.plot(mat['xsall'][iup], mat['cp_steady_data'][iup,2], "o",label = "Exp. up",)
        axx.plot(mat['xsall'][iup], mat['cp_steady_data'][iup,3], "o",label = "Exp. lo",)
        axx.set_xlabel("$x$")
        axx.set_ylabel("$c_p(0)$")
        axx.set_ylim(-1.0, 1.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'{ysmes[i]/bhalf:.3f}')
        plt.show()
    # End for
    
    # Calculate the unsteady pressure coefficients
    cpalpha,cphdot,cpalphadot,cph2dot,cpalpha2dot=SDPMcalcs.unsteadysolve_pitchplunge(body,allbodies,Aphi,Bphi,Cphi,barUinf,barVinf,barWinf,k,c0,Mach,beta,cp_order,xf0,yf0,zf0,install_dir)
    # Calculate total oscillatory pressure distribution
    cptot=(cpalpha+1j*k*cpalphadot+(1j*k)**2*cpalpha2dot)
    # Extract and reshape unsteady pressure on wing only (not wingtips)
    cptot=np.reshape(cptot[0:body['m'][0]*body['n'][0]],(body['m'][0],body['n'][0]),order='C')
    # Calculate pressure jump across the surface
    Deltacptot=np.flipud(cptot[0:body['m'][0]//2,:])-cptot[body['m'][0]//2:body['m'][0],:]
    
    # Interpolate unsteady pressure distribution onto spanwise measurement
    # positions
    cpspan,xspan=SDPMcalcs.cp_interpspan(body,0,cptot,yumes)

    # Plot unsteady pressures
    for i in range (0,len(yumes)):
        iup=np.argwhere(np.round(mat['yuall']*10000.0) == np.round(yumes[i]*10000.0))
        iup=iup[:,0]
        fig, axx = plt.subplots(2,1)
        axx[0].plot(xspan[:,i], np.real(cpspan[:,i]), label = "SDPM")
        axx[0].plot(mat['xuall'][iup], np.real(mat['cp_unsteady_data'][iup,8])/alpha1, "o",label = "Exp. up",)
        axx[0].plot(mat['xuall'][iup], -np.real(mat['cp_unsteady_data'][iup,9])/alpha1, "o",label = "Exp. lo",)
        axx[0].set_xlabel("$x$")
        axx[0].set_ylabel("$\Re(c_p(k))/\\alpha_1$")
        axx[0].set_ylim(-15.0, 5.0)
        axx[0].grid()
        axx[0].set_title("$M_{\\infty}=$" +str(Mach) + ", $\\alpha_0$=" + str(alpha0*180/np.pi) + "$^o$, $2y/b$=" + f'{yumes[i]/bhalf:.3f}')
        axx[0].legend(loc="lower right")
        axx[1].plot(xspan[:,i], np.imag(cpspan[:,i]), label = "SDPM")
        axx[1].plot(mat['xuall'][iup], np.imag(mat['cp_unsteady_data'][iup,8])/alpha1, "o",label = "Exp. up",)
        axx[1].plot(mat['xuall'][iup], -np.imag(mat['cp_unsteady_data'][iup,9])/alpha1, "o",label = "Exp. lo",)
        axx[1].set_xlabel("$x$")
        axx[1].set_ylabel("$\Im(c_p(k))/\\alpha_1$")
        axx[1].set_ylim(-5.0, 5.0)
        axx[1].grid()
        plt.show()
    # End for   
 
# End for
