diff --git a/pyavlpackage/CMakeLists.txt b/pyavlpackage/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ba9acde4f422aa40fde43f4759ff4653b1a29483
--- /dev/null
+++ b/pyavlpackage/CMakeLists.txt
@@ -0,0 +1,6 @@
+# Add the package to the package list for exporting the target
+# and propagate the resulting list back to the parent scope
+list( APPEND PYTHON_TARGETS ${CMAKE_CURRENT_LIST_DIR} )
+set( PYTHON_TARGETS ${PYTHON_TARGETS} PARENT_SCOPE )
+
+install(DIRECTORY ${CMAKE_CURRENT_LIST_DIR} DESTINATION lib)
diff --git a/pyavlpackage/README.md b/pyavlpackage/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..5a581c8d007c8e76063264659deb3eb2c51b4655
--- /dev/null
+++ b/pyavlpackage/README.md
@@ -0,0 +1,11 @@
+# PYAVL Package - Python package for working with AVL (Athena vortex lattice)
+
+## Usage
+Used for Interfacing AVL
+
+## License
+This project is licensed under the GNU General Public License, Version 3
+- see the License.md file for details.
+
+## Contact
+For questions or feedback, please contact contacts@unicado.io
diff --git a/pyavlpackage/pyavl/__init__.py b/pyavlpackage/pyavl/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..eb989e1ef2ad5831b6b59efed2649917e9d38f3b
--- /dev/null
+++ b/pyavlpackage/pyavl/__init__.py
@@ -0,0 +1,24 @@
+# UNICADO - UNIversity Conceptual Aircraft Design and Optimization
+#
+# Copyright (C) 2025 UNICADO consortium
+#
+# 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/>.
+#
+# Description:
+# This file is part of UNICADO.
+
+from .pyavlcore import *
+from .pyavlfileread import *
+from .pyavlfilewrite import *
+from .pyavlrunner import *
\ No newline at end of file
diff --git a/pyavlpackage/pyavl/pyavlcore.py b/pyavlpackage/pyavl/pyavlcore.py
new file mode 100644
index 0000000000000000000000000000000000000000..645ea2cb2b092a4b0ab72553e830f24e8d9b7d00
--- /dev/null
+++ b/pyavlpackage/pyavl/pyavlcore.py
@@ -0,0 +1,423 @@
+# UNICADO - UNIversity Conceptual Aircraft Design and Optimization
+#
+# Copyright (C) 2025 UNICADO consortium
+#
+# 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/>.
+#
+# Description:
+# This file is part of UNICADO.
+
+"""
+PyAvl Core classes
+"""
+from .pyavlfilewrite import *
+from .pyavlfileread import *
+from ambiance import Atmosphere
+
+"""
+Geometry classes (.avl)
+"""
+
+
+class Header:
+ def __init__(self, runcase: str="base", mach: float=0.0, sym: list=[0, 0, 0.0], base: list[float]=[1.0, 1.0, 1.0], ref: list[float]=[0.0, 0.0, 0.0]):
+ """ Initialize avl header - please see avl primer
+
+ Args:
+ runcase (str, optional): Runcase name. Defaults to "base".
+ mach (float, optional): mach number. Defaults to 0.0.
+ sym (list, optional): symmetry elements according to AVL. Defaults to [0, 0, 0.0].
+ base (list[float], optional): Basic aerodynamic parameters (Sref, cref, bref). Defaults to [1.0, 1.0, 1.0].
+ ref (list[float], optional): Reference cg point. Defaults to [0.0, 0.0, 0.0].
+ """
+ self.runcase = runcase
+ self.mach = mach
+ self.sym = sym
+ self.base = base
+ self.ref = ref
+ self.author = ""
+
+ def write(self, fh) -> None:
+ """ Write data to file via fh -> filehandler
+
+ Args:
+ fh (): filehandler/filedescriptor
+ """
+ write_separation_big(fh, "BASE FILE")
+ write_headline(fh, "Geometric Name")
+ write_elem_tab(fh, self.runcase, 2)
+ write_headline(fh, "Mach (normal to wing)")
+ write_elem_tab(fh, self.mach, 2)
+ write_headline(fh, "Symmetric properties (iYsym, iZsym, Zsym)")
+ write_elem_tab(fh, self.sym, 2)
+ write_headline(fh, "Reference properties Sref, Cref, Bref")
+ write_elem_tab(fh, self.base, 2)
+ write_headline(fh, "Reference CG position (X,Y,Z)")
+ write_elem_tab(fh, self.ref, 2)
+ write_separation_big(fh, "Begin of Geometry")
+
+class Control:
+ def __init__(self, tag="ctrl", gain=1.0, xhinge=0.7, xyzhvec=[0.0, 1.0, 0.0], sgndup=1.0):
+ self.key = "CONTROL"
+ self.tag = tag
+ self.gain = gain
+ self.xhinge = xhinge
+ self.xyzhvec = xyzhvec
+ self.sgndup = sgndup
+ self.description = ["!name", "gain", "xhinge", "XYZhvec", "SgnDup"]
+
+ def write(self, fh) -> None:
+ """ Write data to file via fh -> filehandler
+
+ Args:
+ fh (): filehandler/filedescriptor
+ """
+ write_key(fh, self.key)
+ write_elem_tab(fh, self.description, 1)
+ write_elem_tab(fh, self.tag)
+ write_elem_tab(fh, self.gain)
+ write_elem_tab(fh, self.xhinge)
+ write_elem_tab(fh, self.xyzhvec)
+ write_elem_tab(fh, self.sgndup, 2)
+
+
+class Section:
+ def __init__(self, refLE: list[float]=[0.0, 0.0, 0.0], chord: float=1.0, dihedral: float=0.0, nspan: int=1, sspace: float=1.0, afileKey: str="NACA",
+ afile: str="0012"):
+ """ Section class
+
+ Args:
+ refLE (list[float], optional): Position of section reference point (leading edge). Defaults to [0.0, 0.0, 0.0].
+ chord (float, optional): chord length. Defaults to 1.0.
+ dihedral (float, optional): dihedral. Defaults to 0.0.
+ nspan (int, optional): number of spanwise horseshoe vortices. Defaults to 1.
+ sspace (float, optional): spanwise vortex spacing. Defaults to 1.0.
+ afileKey (str, optional): airfoil file key (NACA or AFILE). Defaults to "NACA".
+ afile (str, optional): for key NACA -> enter 4 digit naca profile, otherwise if AFILE -> enter path to airfoil. Defaults to "0012".
+ """
+ self.key = "SECTION"
+ self.ref = refLE
+ self.chord = chord
+ self.dihedral = dihedral
+ self.nspan = nspan
+ self.sspace = sspace
+ self.afileKey = afileKey
+ self.afile = afile
+ self.controls = []
+ self.description = ["!Xle", "Yle", "Zle", "Chord", "dihedral", "nspan", "sspace"]
+
+ def add_control_device(self, control: Control) -> None:
+ """
+
+ Args:
+ control (Control): Add control object
+ """
+ self.controls.append(control)
+
+ def write(self, fh) -> None:
+ """ Write data to file via fh -> filehandler
+
+ Args:
+ fh (): filehandler/filedescriptor
+ """
+ write_separation_small(fh, "Section")
+ write_key(fh, self.key)
+ write_elem_tab(fh, self.description, 1)
+ write_elem_tab(fh, self.ref)
+ write_elem_tab(fh, self.chord)
+ write_elem_tab(fh, self.dihedral)
+ write_elem_tab(fh, self.nspan)
+ write_elem_tab(fh, self.sspace, 2)
+ write_afile(fh, self.afileKey, self.afile)
+ write_new_line(fh)
+ for control in self.controls:
+ control.write(fh)
+
+
+class Yduplicate:
+ def __init__(self, ydupl: float=0.0):
+ """ Yduplicate
+
+ Args:
+ ydupl (float, optional): Y postion of X-Z plane -> only use if iYsym is 0. Defaults to 0.0.
+ """
+ self.key = "YDUPLICATE"
+ self.ydupl = ydupl
+ self.description = ["!YDuplicate - Geometric Symmetry, non aerodynamic symmetry"]
+
+ def write(self, fh) -> None:
+ """ Write yduplicate
+
+ Args:
+ fh (): filehandler/filedescriptor
+ """
+ write_key(fh, self.key)
+ write_elem_tab(fh, self.description, 1)
+ write_elem_tab(fh, self.ydupl, 2)
+
+
+class Surface:
+ def __init__(self, tag: str="surf", nchord: int=10, cspace: float=1.0, nspan: int=20, sspace: float=1.0):
+ """ General surface
+
+ Args:
+ tag (str, optional): Surface tag. Defaults to "surf".
+ nchord (int, optional): Number of chordwise horseshoe vortices. Defaults to 10.
+ cspace (float, optional): Chordwise vortex spacing parameter. Defaults to 1.0.
+ nspan (int, optional): Number of spanwise horseshoe vortices. Defaults to 20.
+ sspace (float, optional): Spanwise vortex spacing. Defaults to 1.0.
+ """
+ self.key = "SURFACE"
+ self.tag = tag
+ self.nchord = nchord
+ self.cspace = cspace
+ self.nspan = nspan
+ self.sspace = sspace
+ self.yduplicates = []
+ self.angles = []
+ self.scales = []
+ self.translates = []
+ self.sections = []
+
+ def add_yduplicate(self, ydupl) -> None:
+ """ Adds yduplicate object
+
+ Args:
+ ydupl (object): yduplicate object
+ """
+ self.yduplicates.append(ydupl)
+
+ def add_angle(self, angle) -> None:
+ """ Adds angle object
+
+ Args:
+ angle (object): angle object
+ """
+ self.angles.append(angle)
+
+ def add_scale(self, scale) -> None:
+ """ Adds scale object
+
+ Args:
+ scale (object): scale object
+ """
+ self.scales.append(scale)
+
+ def add_translate(self, translate) -> None:
+ """ Add translate object
+
+ Args:
+ translate (object): translate object
+ """
+ self.translates.append(translate)
+
+ def add_section(self, section) -> None:
+ """ Add section object
+
+ Args:
+ section (object): section object
+ """
+ self.sections.append(section)
+
+ def write(self, fh) -> None:
+ """ Write data to file via fh -> filehandler
+
+ Args:
+ fh (): filehandler/filedescriptor
+ """
+ write_separation_big(fh, f'Begin of Surface [{self.tag}]')
+ write_key(fh, self.key)
+ write_tag(fh, self.tag)
+ write_elem_tab(fh, self.nchord)
+ write_elem_tab(fh, self.cspace)
+ write_elem_tab(fh, self.nspan)
+ write_elem_tab(fh, self.sspace, 2)
+
+ for ydl in self.yduplicates:
+ ydl.write(fh)
+ for angle in self.angles:
+ angle.write(fh)
+ for scale in self.scales:
+ scale.write(fh)
+ for translate in self.translates:
+ translate.write(fh)
+ for section in self.sections:
+ section.write(fh)
+
+ write_separation_big(fh, f'End of Surface [{self.tag}]')
+
+
+class Scale:
+ def __init__(self, xscale: float=1.0, yscale: float=1.0, zscale: float=1.0):
+ """ Scale object
+
+ Args:
+ xscale (float, optional): x-scale factor. Defaults to 1.0.
+ yscale (float, optional): y-scale factor. Defaults to 1.0.
+ zscale (float, optional): z-scale factor. Defaults to 1.0.
+ """
+ self.key = "SCALE"
+ self.xscale = xscale
+ self.yscale = yscale
+ self.zscale = zscale
+ self.description = ["!Xscale", "Yscale", "Zscale"]
+
+ def write(self, fh) -> None:
+ """ Write to file
+
+ Args:
+ fh (): filehandler/filedescriptor
+ """
+ write_key(fh, self.key)
+ write_elem_tab(fh, self.description, 1)
+ write_elem_tab(fh, self.xscale)
+ write_elem_tab(fh, self.yscale)
+ write_elem_tab(fh, self.zscale, 2)
+
+
+class Translate:
+ def __init__(self, dX: float=0.0, dY: float=0.0, dZ: float=0.0):
+ """ Translate object (right hand coordinate system starting in nose through back (x),
+ left (y), up (z))
+
+ Args:
+ dX (float, optional): x translation. Defaults to 0.0.
+ dY (float, optional): y translation. Defaults to 0.0.
+ dZ (float, optional): z translation. Defaults to 0.0.
+ """
+ self.key = "TRANSLATE"
+ self.dX = dX
+ self.dY = dY
+ self.dZ = dZ
+ self.description = ["!dX", "dY", "dZ"]
+
+ def write(self, fh) -> None:
+ """ Write to file
+
+ Args:
+ fh (): filehandler/filedescriptor
+ """
+ write_key(fh, self.key)
+ write_elem_tab(fh, self.description, 1)
+ write_elem_tab(fh, self.dX)
+ write_elem_tab(fh, self.dY)
+ write_elem_tab(fh, self.dZ, 2)
+
+
+class Angle:
+ def __init__(self, dAinc: float=0.0) -> None:
+ """ Angle object
+
+ Args:
+ dAinc (float, optional): Offset added on to the Ainc values (deg). Defaults to 0.0.
+ """
+ self.key = "ANGLE"
+ self.dAinc = dAinc
+ self.description = ["!Dihedral"]
+
+ def write(self, fh) -> None:
+ """ Write to file
+
+ Args:
+ fh (): filehandler/filedescriptor
+ """
+ write_key(fh, self.key)
+ write_elem_tab(fh, self.description, 1)
+ write_elem_tab(fh, self.dAinc, 1)
+
+
+"""
+Runcase class (.run)
+"""
+
+
+class Runcase:
+ def __init__(self, alpha: list[float]=[0.0], beta: float=0.0, mach: float=0.0, cgID: int=1, cgpos: list[float]=[0.0, 0.0, 0.0], height: float=11000.0,
+ aileron: float=0.0,
+ elevator: float=0.0, rudder: float=0.0):
+ """ Initialize a single runcase for avl -> inputs for .run file
+
+ Args:
+ alpha (list[float], optional): angle of attack. Defaults to [0.0].
+ mach (float, optional): mach number. Defaults to 0.0.
+ beta (float, optional): angle. Defaults to 0.0.
+ cgID (int, optional): center of gravity id (depends on user). Defaults to 1.
+ cgpos (list[float], optional): cgpos w.r.t. cgID. Defaults to [0.0, 0.0, 0.0].
+ height (list[float], optional): height/altitude. Defaults to [11000.0].
+ aileron (float, optional): aileron (roll) control surface deflection. Defaults to 0.0.
+ elevator (float, optional): elevator (pitch) control surface deflection. Defaults to 0.0.
+ rudder (float, optional): rudder (yaw) control surface deflection. Defaults to 0.0.
+ """
+ atmo = Atmosphere(height)
+
+ if alpha is None:
+ alpha = [0.0]
+ self.alpha = alpha
+ self.beta = beta
+ self.mach = mach
+ self.aileron = aileron
+ self.elevator = elevator
+ self.rudder = rudder
+ self.cgID = cgID
+ self.cgpos = cgpos
+ self.density = float(atmo.density[0])
+ self.a = float(atmo.speed_of_sound[0])
+ self.numOfRuncases = len(alpha)
+ self.cases = []
+ self.caseName = None
+
+ def create_runcase_fileName(self) -> str:
+ """ Create runcase filename
+
+ Returns:
+ str -- runcase filename
+ """
+ strmach = f'mach{self.mach}_'
+ strcg = f'cgID{self.cgID}'
+ s = "Runcase_" + strmach + strcg
+ s = s.split('.')
+ s = '_'.join(s)
+ return s
+
+ def create_runcase_file(self,id: int) -> None:
+ """ Creates a .run file based on a Runcase
+
+ Args:
+ id (int): Runcase file id
+ """
+ self.caseName = self.create_runcase_fileName()
+ idx = 1
+ with open(self.caseName + ".run", "w") as file:
+ for a in self.alpha:
+ self.cases.append(f'file_id{id+idx}')
+ write_runcase(file, idx, self.cases[-1], alpha=a, beta=self.beta, mach=self.mach, density=self.density,
+ a=self.a, aileron=self.aileron, elevator=self.elevator,
+ rudder=self.rudder, cgref=self.cgpos)
+ idx += 1
+
+ def get_runcases(self) -> list:
+ """ Returns list of runcases
+
+ Returns:
+ list: list of runcases
+ """
+ return self.cases
+
+ def get_runcase_filename(self) -> str:
+ """ Get runcase filename
+
+ Returns:
+ str: runcase filename
+ """
+ return self.caseName
diff --git a/pyavlpackage/pyavl/pyavlfileread.py b/pyavlpackage/pyavl/pyavlfileread.py
new file mode 100644
index 0000000000000000000000000000000000000000..48a6fa0f25a178c4f63c816d05006f05ddbfca24
--- /dev/null
+++ b/pyavlpackage/pyavl/pyavlfileread.py
@@ -0,0 +1,321 @@
+# UNICADO - UNIversity Conceptual Aircraft Design and Optimization
+#
+# Copyright (C) 2025 UNICADO consortium
+#
+# 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/>.
+#
+# Description:
+# This file is part of UNICADO.
+
+"""
+Python AVL Class for Reading Data
+
+This module provides classes to process AVL stability files and extract relevant aerodynamics data.
+"""
+import numpy as np
+import re
+import csv
+import math
+
+class PyAvlStabilityFilesConvert:
+ """
+ Converts AVL stability files to a structured format and writes them to a CSV file.
+
+ Attributes:
+ filenames (list): List of AVL stability filenames.
+ datatable (numpy.ndarray): Table of extracted data.
+ keys (list): List of keys corresponding to the data columns.
+ stabfiles (list): List of PyAvlStabilityFile objects.
+ keyListControls (list, optional): List of control keys.
+ """
+ def __init__(self, filenames=[], keyListControls=None, fileout="default.csv"):
+ """
+ Initializes the PyAvlStabilityFilesConvert class and processes the given stability files.
+
+ Args:
+ filenames (list): List of filenames to process.
+ keyListControls (list, optional): List of control keys.
+ fileout (str): Output CSV filename.
+ """
+ filenames.sort()
+ self.filenames = filenames
+ self.datatable = None
+ self.keys = None
+ self.stabfiles = []
+ self.keyListControls = keyListControls
+ for filename in self.filenames:
+ self.stabfiles.append(PyAvlStabilityFile(filename, keylistControls=keyListControls))
+
+ with open(fileout, "w", newline='') as file:
+ if self.keys is None:
+ self.keys = self.stabfiles[0].get_keys()
+
+ writer = csv.DictWriter(file, fieldnames=self.keys)
+ writer.writeheader()
+ datatable = []
+ for stabfile in self.stabfiles:
+ writer.writerow(stabfile.get_data())
+ tmp = stabfile.get_data()
+ datatable.append(list(tmp.values()))
+ self.datatable = np.array(datatable)
+
+ def get_keys(self):
+ """Returns the list of keys from the stability files."""
+ return self.keys
+
+ def get_data(self):
+ """Returns the numerical data extracted from the stability files."""
+ return self.datatable
+
+ def get_data_by_key(self,key):
+ """
+ Retrieves data for a specific key.
+
+ Args:
+ key (str): Key for which data is requested.
+
+ Returns:
+ numpy.ndarray or None: Data for the given key, or None if not found.
+ """
+ try:
+ key_index = self.keys.index(key)
+ return self.datatable[:,key_index]
+ except ValueError:
+ return None
+
+
+ def get_data_by_control_key(self,key: str ,controlkey: str, exact: bool=False):
+ """
+ Retrieves data based on control key modifications.
+
+ Args:
+ key (str): The aerodynamic coefficient key.
+ controlkey (str): The control surface key.
+ exact (bool, optional): Whether to look for an exact match of the control key.
+
+ Returns:
+ numpy.ndarray or None: The extracted data.
+ """
+ available_keys = ["CL", "CY", "Cl", "Cm", "Cn","CDff"]
+ try:
+ # Get key in available key list -> set to 1 since zero would be direct control key
+ available_key_index = available_keys.index(key) + 1
+ if not exact:
+ control_key_indices = [idx + available_key_index for idx, item in enumerate(self.keys) if item.startswith(controlkey)]
+ # if more than one key index is available
+ if len(control_key_indices) > 1:
+ # Combine data
+ return np.sum([self.datatable[:, idx]*math.degrees(1) for idx in control_key_indices], axis=1)
+ else:
+ # else
+ return self.datatable[:, control_key_indices[0]]*math.degrees(1)
+ else:
+ control_key_index = self.keys.index(controlkey) + available_key_index
+ return self.datatable[:, control_key_index]*math.degrees(1)
+
+
+ except ValueError:
+ return None
+
+
+
+
+class PyAvlStabilityFiles:
+ """
+ Manages multiple AVL stability files.
+
+ Attributes:
+ filenames (list): List of AVL stability filenames.
+ stabfiles (list): List of processed stability file data.
+ keys (list): List of data keys.
+ ctrlKeys (list, optional): List of control keys.
+ """
+ def __init__(self,filenames=[], keyListControls=None):
+ """
+ Initializes the class and processes the given stability files.
+
+ Args:
+ filenames (list): List of filenames.
+ keyListControls (list, optional): List of control keys.
+ """
+ self.filenames = filenames
+
+ self.stabfiles = []
+ self.keys = None
+ self.ctrlKeys = keyListControls
+ for filename in self.filenames:
+ stabilityfile = PyAvlStabilityFile(filename, keylistControls=keyListControls)
+ if self.keys is None:
+ self.keys = stabilityfile.get_keys()
+ self.stabfiles.append(stabilityfile.get_data())
+
+ def get_data(self):
+ """Returns the data extracted from the stability files."""
+ return self.stabfiles
+
+ def get_keys(self):
+ """Returns the list of keys from the stability files."""
+ return self.keys
+
+class PyAvlStabilityFile:
+ """
+ Represents an AVL stability file.
+
+ Attributes:
+ filename (str): Name of the AVL stability file.
+ keylist (list): List of extracted data keys.
+ keylistControls (list): List of control keys.
+ data (dict): Extracted data from the file.
+ """
+
+ def __init__(self, filename, keylistControls=None):
+ """
+ Initializes the stability file by parsing its contents.
+
+ Args:
+ filename (str): Filename of the AVL stability file.
+ keylistControls (list, optional): List of control keys.
+ """
+ self.filename = filename
+
+ self.keylist = self.init_dictionary()
+ self.keylistControls = self.set_control_keys(keylistControls)
+ self.keyControlNums = []
+ self.data = self.fill_dictionary()
+
+ def get_keys(self):
+ """Returns the list of keys in the stability file."""
+ return self.keylist
+
+ def get_data(self):
+ """Returns the extracted data from the stability file."""
+ return self.data
+
+ def init_dictionary(self):
+ """Initializes the dictionary with predefined aerodynamic data keys."""
+ keylist = []
+
+ keylist += ["Altitude","Sref", "Cref", "Bref"]
+ keylist += ["Xref", "Yref", "Zref", "Xnp"]
+ keylist += ["Alpha", "Beta", "Mach", "pb/2V", "qc/2V", "rb/2V"]
+ keylist += ["CLtot","CYtot", "CDtot", "CDvis", "CDind", "Cltot", "Cmtot", "Cntot"]
+
+ namederiv = ["CL", "CY", "Cl", "Cm", "Cn"]
+ addderiv = ["a", "b", "p", "q", "r"]
+
+ for id in addderiv:
+ keylist += [val + id for val in namederiv]
+ return keylist
+
+ def fill_dictionary(self):
+ """
+ Reads the AVL stability file and fills the data dictionary.
+
+ Returns:
+ dict: Extracted aerodynamic data mapped to their respective keys.
+ """
+ with open(self.filename, "r") as file:
+ lines = file.readlines()
+
+ addderiv = []
+ namederiv = ["CL", "CY", "Cl", "Cm", "Cn","CDff"]
+ for ctrlkey in self.keylistControls:
+ id = self.get_control_device_num(lines, ctrlkey)
+ self.keylist.append(ctrlkey)
+ self.keylist += [val + id for val in namederiv]
+
+ data = dict.fromkeys(self.keylist, None)
+
+ for key in self.keylist:
+ if key == "Altitude":
+ data["Altitude"] = self.get_ft_from_runcase(lines)
+ else:
+ data[key] = self.get_x_value(lines, key)
+
+ return data
+
+ def set_control_keys(self, controlnames=None):
+ """Sets the control keys."""
+ controlnames = list(controlnames)
+ keylistControls = []
+ if controlnames is None or type(controlnames) is not type([]):
+ return keylistControls
+ else:
+ keylistControls = controlnames
+ return keylistControls
+
+ def get_x_value(self, lines, key):
+ """
+ Extracts the numerical value associated with a given key from the file.
+
+ Args:
+ lines (list): Lines from the AVL stability file.
+ key (str): The key whose value needs to be extracted.
+
+ Returns:
+ float or None: Extracted numerical value or None if not found.
+ """
+ searchObj = re.compile(rf'{key} * = *.\d*.\d*')
+ for line in lines:
+ res = searchObj.search(line)
+ if res is not None:
+ ret = res.group().split("=")
+ return float(ret[1])
+
+ def get_control_device_num(self, lines, ctrlkey):
+ """
+ Retrieves the control device number for a given control key.
+
+ Args:
+ lines (list): Lines from the AVL stability file.
+ ctrlkey (str): The control key whose number is needed.
+
+ Returns:
+ str: Control device number as a string.
+ """
+ searchObj = re.compile(rf'{ctrlkey[:]} * d\d*.\d*')
+ for line in lines:
+ res = searchObj.search(line)
+ if res is not None:
+ ret = res.group().split()
+ return ret[1]
+
+ def get_ft_from_runcase(self,lines):
+ """
+ Extracts altitude information from the run case.
+
+ Args:
+ lines (list): Lines from the AVL stability file.
+
+ Returns:
+ float: Altitude value extracted from the run case.
+ """
+ searchObj = re.compile(rf'Run case: RunCase \d* - \d*.\d*')
+ for line in lines:
+ res = searchObj.search(line)
+ if res is not None:
+ ret = res.group().split("-")
+ return float(ret[-1])
+
+ def write(self, fileout):
+ """
+ Writes the extracted data to a CSV file.
+
+ Args:
+ fileout (str): The output filename without extension.
+ """
+ with open(fileout + ".csv", "w", newline='') as file:
+ writer = csv.DictWriter(file, fieldnames=self.keylist)
+ writer.writeheader()
+ writer.writerow(self.data)
diff --git a/pyavlpackage/pyavl/pyavlfilewrite.py b/pyavlpackage/pyavl/pyavlfilewrite.py
new file mode 100644
index 0000000000000000000000000000000000000000..45def28a7c525a17017db9af2545918e05a15e7a
--- /dev/null
+++ b/pyavlpackage/pyavl/pyavlfilewrite.py
@@ -0,0 +1,233 @@
+# UNICADO - UNIversity Conceptual Aircraft Design and Optimization
+#
+# Copyright (C) 2025 UNICADO consortium
+#
+# 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/>.
+#
+# Description:
+# This file is part of UNICADO.
+
+"""
+Python AVL Class for Reading Data
+
+This module provides classes to process AVL stability files and extract relevant aerodynamics data.
+"""
+
+def write_elem_new_line(fh, val):
+ """
+ Writes a value to the file and adds a new line.
+
+ Args:
+ fh (file object): File handle.
+ val (str): Value to write.
+ """
+ if len(val):
+ fh.write(f'{val}\n')
+
+
+def write_elem_tab(fh, elems, newline=0):
+ """
+ Writes elements separated by spaces and optionally adds a new line.
+
+ Args:
+ fh (file object): File handle.
+ elems (list or str): Elements to write.
+ newline (int, optional): Number of new lines to add.
+ """
+ if type(elems) is list:
+ for elem in elems:
+ fh.write(f'{elem} ')
+ else:
+ fh.write(f'{elems} ')
+
+ write_new_line(fh, newline)
+
+
+def write_key(fh, key):
+ """Writes a key as a new line."""
+ write_elem_new_line(fh, key)
+
+
+def write_tag(fh, tag):
+ """Writes a tag as a new line."""
+ write_elem_new_line(fh, tag)
+
+
+def write_afile(fh, key, afile):
+ """
+ Writes a key followed by a filename as new lines.
+
+ Args:
+ fh (file object): File handle.
+ key (str): Key to write.
+ afile (str): Associated filename.
+ """
+ write_elem_new_line(fh, key)
+ write_elem_new_line(fh, f'{afile}')
+
+
+def write_commentline(list, commentstring):
+ """
+ Adds a commented line to a list.
+
+ Args:
+ lst (list): List to append the comment.
+ commentstring (str): Comment text.
+ """
+ commentline = "# " + commentstring + "\n"
+ list.append(commentline)
+
+
+def write_separation_big(file, token=""):
+ """
+ Writes a large separation line with an optional token.
+
+ Args:
+ file (file object): File handle.
+ token (str, optional): Text to include in the separation.
+ """
+ sepString = "# =="
+ if len(token):
+ sepString += f' {token} '
+
+ lsepString = len(sepString)
+ for i in range(0, 80 - lsepString):
+ sepString += "="
+ sepString += "\n"
+ file.write(sepString)
+
+
+def write_separation_small(file, token=""):
+ """
+ Writes a small separation line with an optional token.
+
+ Args:
+ file (file object): File handle.
+ token (str, optional): Text to include in the separation.
+ """
+ sepString = "# --"
+ if len(token):
+ sepString += f' {token} '
+
+ lsepString = len(sepString)
+ for i in range(0, 80 - lsepString):
+ sepString += "-"
+ sepString += "\n"
+ file.write(sepString)
+
+
+def write_headline(file, short=""):
+ """Writes a headline in a formatted manner."""
+ file.write(f'# -- {short} --\n')
+
+
+def write_new_line(fh, num=1):
+ """
+ Writes the specified number of new lines.
+
+ Args:
+ fh (file object): File handle.
+ num (int): Number of new lines to write.
+ """
+ crnl = ""
+ for i in range(0, num):
+ crnl += "\n"
+ fh.write(crnl)
+
+
+def write_avl_file(author="", datalist=[], filename="filename", postfix=".avl"):
+ """
+ Writes AVL data to a specified file.
+
+ Args:
+ author (str, optional): Author name.
+ datalist (list): Data elements to write.
+ filename (str): Base filename.
+ postfix (str, optional): File extension.
+ """
+ with open(filename + postfix, "w") as avlFile:
+ for item in datalist:
+ item.write(avlFile)
+
+
+def write_runcase(file, runcasenum=1, runcasename="default", alpha=0.0, beta=0.0, pb2v=0.0, qc2v=0.0, rb2v=0.0,
+ controls=None, mach=0.0, CD0=0.02, density=1.225,a=340.294, cgref=[0.0, 0.0, 0.0]):
+ """
+ Writes a run case block to the file.
+
+ Args:
+ file (file object): File handle.
+ runcasenum (int, optional): Run case number.
+ runcasename (str, optional): Run case name.
+ alpha (float, optional): Angle of attack.
+ beta (float, optional): Sideslip angle.
+ pb2v, qc2v, rb2v (float, optional): Non-dimensional rotational rates.
+ controls (list, optional): List of control surfaces.
+ mach (float, optional): Mach number.
+ CD0 (float, optional): Zero-lift drag coefficient.
+ density (float, optional): Air density.
+ a (float, optional): Speed of sound.
+ cgref (list, optional): Center of gravity reference coordinates.
+ """
+ file.write(f' ---------------------------------------------\n')
+ file.write(f' Run case {runcasenum}: {runcasename}\n')
+ file.write(f' alpha -> alpha = {alpha}\n')
+ file.write(f' beta -> beta = {beta}\n')
+ file.write(f' pb/2V -> pb/2V = {pb2v}\n')
+ file.write(f' qc/2V -> qc/2V = {qc2v}\n')
+ file.write(f' rb/2V -> rb/2V = {rb2v}\n')
+ if controls is not None and isinstance(controls,list):
+ for control in controls:
+ file.write(f' {control}\t\t->{control}\t\t=\t{0.0}\n')
+ # if aileron is not None:
+ # file.write(f' aileron -> aileron = {aileron}\n')
+ # if elevator is not None:
+ # file.write(f' elevator -> elevator = {elevator}\n')
+ # if rudder is not None:
+ # file.write(f' rudder -> rudder = {rudder}\n')
+
+ file.write("\n")
+
+ default = 0.0
+ velocity = mach*a[0]
+ file.write(f' alpha = {alpha} deg\n')
+ file.write(f' beta = {beta} deg\n')
+ file.write(f' pb/2V = {default}\n')
+ file.write(f' qc/2V = {default}\n')
+ file.write(f' rb/2V = {default}\n')
+ file.write(f' CL = {default}\n')
+ file.write(f' CDo = {CD0}\n')
+ file.write(f' bank = {default} deg\n')
+ file.write(f' elevation = {default} deg\n')
+ file.write(f' heading = {default} deg\n')
+ file.write(f' Mach = {mach}\n')
+ file.write(f' velocity = {velocity:.4f} Lunit/Tunit\n')
+ file.write(f' density = {density[0]:.4f} Munit/Lunit^3\n')
+ file.write(f' grav.acc. = {9.80665} Lunit/Tunit^2\n')
+ file.write(f' turn_rad. = {default} Lunit\n')
+ file.write(f' load_fac. = {1.00000}\n')
+ file.write(f' X_cg = {cgref[0]} Lunit\n')
+ file.write(f' Y_cg = {cgref[1]} Lunit\n')
+ file.write(f' Z_cg = {cgref[2]} Lunit\n')
+ file.write(f' mass = {1.00000} Munit\n')
+ file.write(f' Ixx = {1.00000} Munit-Lunit^2\n')
+ file.write(f' Iyy = {1.00000} Munit-Lunit^2\n')
+ file.write(f' Izz = {1.00000} Munit-Lunit^2\n')
+ file.write(f' Ixy = {0.00000} Munit-Lunit^2\n')
+ file.write(f' Iyz = {0.00000} Munit-Lunit^2\n')
+ file.write(f' Izx = {0.00000} Munit-Lunit^2\n')
+ file.write(f' visc CL_a = {0.00000}\n')
+ file.write(f' visc CL_u = {0.00000}\n')
+ file.write(f' visc CM_a = {0.00000}\n')
+ file.write(f' visc CM_u = {0.00000}\n\n')
diff --git a/pyavlpackage/pyavl/pyavlrunner.py b/pyavlpackage/pyavl/pyavlrunner.py
new file mode 100644
index 0000000000000000000000000000000000000000..a2cb5a2a346168d94465a1c616a0ae58d3d19b25
--- /dev/null
+++ b/pyavlpackage/pyavl/pyavlrunner.py
@@ -0,0 +1,211 @@
+# UNICADO - UNIversity Conceptual Aircraft Design and Optimization
+#
+# Copyright (C) 2025 UNICADO consortium
+#
+# 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/>.
+#
+# Description:
+# This file is part of UNICADO.
+
+"""
+PyAvl - A Python to AVL conversion script for generating lifting surfaces based on given geometry.
+
+This module provides utilities to run Athena Vortex Lattice (AVL) simulations and install AVL/XFOIL tools.
+"""
+import os
+import subprocess
+import platform
+import urllib.request
+import shutil
+import stat
+
+class PyAvlRunner:
+ """
+ Handles execution of AVL simulations, manages run cases, and processes stability results.
+ """
+ def __init__(self,avlFile=None,runcasefile=None,runcases=["default"]):
+ """
+ Initializes the PyAvlRunner class.
+
+ Args:
+ avlFile (str, optional): Path to AVL geometry file.
+ runcasefile (str, optional): Path to AVL run case file.
+ runcases (list, optional): List of run case names.
+ """
+ self.avlproc = None
+ self.cmdline = ""
+ self.avlFile = avlFile
+ self.runcaseFile = runcasefile
+ self.runcases = runcases
+ self.count = 1
+ self.path_to_session_folder = None
+ self.path_to_stability_files = None
+ self.path_to_runcase_files = None
+
+ def session_setup(self):
+ """
+ Sets up the AVL session with geometry and run cases.
+ """
+ # Load geometry and case
+ self.cmdline += f"load {self.avlFile}\n"
+ self.cmdline += f"case {self.runcaseFile}.run\n"
+
+
+ # Switch to operation mode
+ self.cmdline += f"oper\n"
+
+ # Select runcase
+ for id, runcase in enumerate(self.runcases):
+ self.cmdline += f"{id+1}\n"
+ self.cmdline += "x\n"
+ self.cmdline += f"st\n{runcase}_{id}\n"
+ self.cmdline += "\nquit\n"
+
+ def session_setup_single_rc(self, path_to_session_folder, runcase_files_folder, stability_files_folder, results_folder, runcase_ids):
+
+ # Setup paths to session and resulting files
+ self.path_to_session_folder = path_to_session_folder
+ self.path_to_runcase_files = runcase_files_folder
+ self.path_to_stability_files = stability_files_folder
+ self.path_to_session_result_file = results_folder
+
+ # Load geometry and case
+ self.cmdline += f"load {path_to_session_folder}/{self.avlFile}\n"
+ for _, id in enumerate(runcase_ids):
+ self.cmdline += f"case {runcase_files_folder}/rc_{id:03}.run\n"
+
+ # Switch to operation mode
+ self.cmdline += f"oper\n"
+ self.cmdline += f"I\n"
+ self.cmdline += "x\n"
+ self.cmdline += f"st {stability_files_folder}/res_{id:03}.stab\n"
+ self.cmdline += "\n"
+ self.cmdline += "quit\n"
+
+ return self.cmdline
+
+ def run_avl_session(self, cmdline=None):
+ """
+ Executes the AVL session.
+
+ Args:
+ cmdline (str, optional): Custom command sequence for AVL.
+ """
+ if cmdline is None:
+ cmdline = self.cmdline
+ print("Computing aerodynamics with Athena Vortex Lattice (AVL) ... ",end="")
+ process = self.avl_process()
+ process.communicate(input=cmdline.encode())
+ process.wait()
+ print("finished")
+
+
+ def avl_process(self, path_to_avl_executable: str="."):
+ """
+ Initializes the AVL process.
+
+ Args:
+ path_to_avl_executable (str, optional): Path to the AVL executable.
+
+ Returns:
+ subprocess.Popen: AVL process.
+ """
+ avl_executable = "avl"
+ if platform.system() == "Windows":
+ avl_executable = "avl.exe"
+ return subprocess.Popen(args=[f'{path_to_avl_executable}/{avl_executable}'],
+ stdin=subprocess.PIPE,
+ stdout=subprocess.DEVNULL,
+ stderr=subprocess.DEVNULL,
+ bufsize=0)
+
+ def session_End(self):
+ """End avl session"""
+ self.avlproc.terminate()
+ self.avlproc.kill()
+
+ def session_folder(self):
+ """Return session folder path"""
+ return self.path_to_session_folder
+
+ def stability_files_folder(self):
+ """Returns stability files folder path"""
+ return self.path_to_stability_files
+
+ def runcase_files_folder(self):
+ """Returns runcases files folder path"""
+ return self.path_to_runcase_files
+
+ def result_files_folder(self):
+ """Returns result file folder path"""
+ return self.path_to_session_result_file
+
+
+class AVLinstaller:
+ """
+ Downloads and installs AVL based on the user's operating system.
+ """
+ def __init__(self):
+ """Initializes AVLinstaller with appropriate download URLs."""
+ self.avl_url = {"root": "http://web.mit.edu/drela/Public/web/avl/",
+ "Windows": "avl3.40_execs/WIN64/avl.exe",
+ "Darwin-x86_64": "avl3.40_execs/DARWIN64/avl",
+ "Darwin-arm64": "avl3.40_execs/DARWINM1/avl",
+ "Linux": "avl3.40_execs/LINUX64/avl"}
+ self.operating_system: str=self.check_os()
+ self.machine_architecture: str=self.check_architecture()
+ self.avl_download_link: str=None
+
+ def check_os(self) -> str:
+ """Returns the operating system name."""
+ return platform.system()
+
+ def check_architecture(self) -> str:
+ """Returns the system architecture."""
+ return platform.machine()
+
+ def set_avl_download_link(self) -> None:
+ """Determines the correct AVL download link based on the OS and architecture."""
+ if self.operating_system == "Darwin":
+ self.operating_system += "-" + self.machine_architecture
+
+ self.avl_download_link = self.avl_url["root"] + self.avl_url[self.operating_system]
+
+ def is_command_installed(self, command):
+ """Determines if command is installed correctly"""
+ return shutil.which(command) is not None
+
+ def download_avl(self) -> bool:
+ """Downloads and installs AVL if not already installed."""
+ print("Checking AVL ...")
+ executable = "avl"
+ if self.operating_system == "Windows":
+ executable += ".exe"
+ path_to_executable = "./" + executable
+
+ if self.is_command_installed(path_to_executable) or os.path.exists(path_to_executable):
+ print("AVL is already installed.")
+ return True
+ self.set_avl_download_link()
+ print("Downloading AVL ...")
+ print(self.avl_download_link)
+ print(os.getcwd())
+ urllib.request.urlretrieve(self.avl_download_link, os.getcwd()+"/" +executable)
+ if self.is_command_installed(path_to_executable) or os.path.exists(path_to_executable):
+ print("AVL is installed.")
+ print("Changing rights")
+ os.chmod(path_to_executable, stat.S_IRUSR | stat.S_IWUSR | stat.S_IXUSR | stat.S_IRGRP | stat.S_IXGRP | stat.S_IROTH | stat.S_IXOTH)
+ return True
+ else:
+ return False
diff --git a/pyavlpackage/pyproject.toml b/pyavlpackage/pyproject.toml
new file mode 100644
index 0000000000000000000000000000000000000000..9098306081ae1d03ee075359ba77f1bd2f6d7ae4
--- /dev/null
+++ b/pyavlpackage/pyproject.toml
@@ -0,0 +1,27 @@
+[build-system]
+requires = ["setuptools>=42", "wheel"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "pyavl"
+version = "1.0.0"
+description = "Simple python avl interface to build .avl files and runcase files."
+authors = [{ name = "Christopher Ruwisch", email = "christopher.ruwisch@tu-berlin.de" }]
+readme = "README.md"
+license = { file = "LICENSE" }
+keywords = ["avl"]
+classifiers = [
+ "Programming Language :: Python :: 3",
+ "License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)",
+ "Operating System :: OS Independent",
+]
+dependencies = [
+ "numpy>=1.26.4",
+]
+
+[tool.setuptools]
+packages = ["pyavl"]
+
+[project.urls]
+homepage = "https://unicado.pages.rwth-aachen.de/unicado.gitlab.io/"
+repository = "https://git.rwth-aachen.de/unicado/libraries"
\ No newline at end of file
diff --git a/pyavlunicadopackage/CMakeLists.txt b/pyavlunicadopackage/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ba9acde4f422aa40fde43f4759ff4653b1a29483
--- /dev/null
+++ b/pyavlunicadopackage/CMakeLists.txt
@@ -0,0 +1,6 @@
+# Add the package to the package list for exporting the target
+# and propagate the resulting list back to the parent scope
+list( APPEND PYTHON_TARGETS ${CMAKE_CURRENT_LIST_DIR} )
+set( PYTHON_TARGETS ${PYTHON_TARGETS} PARENT_SCOPE )
+
+install(DIRECTORY ${CMAKE_CURRENT_LIST_DIR} DESTINATION lib)
diff --git a/pyavlunicadopackage/README.md b/pyavlunicadopackage/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..965890db22ebc55907404c960fd4ab1f5882bf39
--- /dev/null
+++ b/pyavlunicadopackage/README.md
@@ -0,0 +1,11 @@
+# PYAVLUNICADO Package - Python package for working with AVL (Athena vortex lattice) in combination with Unicado
+
+## Usage
+Used for Interfacing AVL and UNICADO to compute aerodynamic derivatives
+
+## License
+This project is licensed under the GNU General Public License, Version 3
+- see the License.md file for details.
+
+## Contact
+For questions or feedback, please contact contacts@unicado.io
diff --git a/pyavlunicadopackage/pyavlunicado/__init__.py b/pyavlunicadopackage/pyavlunicado/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..fcd206fcb4468e1fd47109ad8a5e5e60b0a8f8cb
--- /dev/null
+++ b/pyavlunicadopackage/pyavlunicado/__init__.py
@@ -0,0 +1,20 @@
+# UNICADO - UNIversity Conceptual Aircraft Design and Optimization
+#
+# Copyright (C) 2025 UNICADO consortium
+#
+# 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/>.
+#
+# Description:
+# This file is part of UNICADO.
+
diff --git a/pyavlunicadopackage/pyavlunicado/avl b/pyavlunicadopackage/pyavlunicado/avl
new file mode 100755
index 0000000000000000000000000000000000000000..01cdfec90ad2f14f496fc4c613d3fcdf489879ac
Binary files /dev/null and b/pyavlunicadopackage/pyavlunicado/avl differ
diff --git a/pyavlunicadopackage/pyavlunicado/avlunicadogeometry.py b/pyavlunicadopackage/pyavlunicado/avlunicadogeometry.py
new file mode 100644
index 0000000000000000000000000000000000000000..b9fc6d1ff6ed3bcb00c4c478115c7abb542e7ee8
--- /dev/null
+++ b/pyavlunicadopackage/pyavlunicado/avlunicadogeometry.py
@@ -0,0 +1,488 @@
+# UNICADO - UNIversity Conceptual Aircraft Design and Optimization
+#
+# Copyright (C) 2025 UNICADO consortium
+#
+# 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/>.
+#
+# Description:
+# This file is part of UNICADO.
+
+"""
+PyAvlUnicado - Geometry Processing for AVL
+
+This module defines classes for handling aerodynamic surfaces and profile scaling within the AVL simulation framework.
+
+Features:
+- Constructs AVL-compatible aerodynamic surface representations.
+- Computes 3D transformations and rotations.
+- Integrates with `pyaircraftgeometry2` for wing and control surface geometry.
+- Scales airfoil profiles for AVL using direct geometric transformations based on XFOIL code.
+
+Dependencies:
+- numpy
+- pyaircraftgeometry2
+- pyavl
+"""
+import pyaircraftgeometry2 as geom2
+import pyavl as avl
+import os
+
+import shutil
+import numpy as np
+
+from collections import defaultdict
+from .xfoilgeometryscale import *
+
+class AerodynamicSurface:
+ """
+ Represents an aerodynamic surface for AVL simulations.
+
+ This class constructs an AVL-compatible aerodynamic surface from aircraft geometry data,
+ allowing for transformations, profile scaling, and control device integration.
+ """
+ def __init__(self, paths, folder_to_create_data: str = ".") -> None:
+ """
+ Initializes an aerodynamic surface object.
+
+ Args:
+ paths (dict): Dictionary containing paths to required data.
+ folder_to_create_data (str, optional): Path where processed airfoil data will be stored. Defaults to "./".
+ """
+ self.aircraft_name = None
+ self.surface = None
+ self.paths = paths
+ self.rotation_geom2_to_avl = None
+ self.folder_to_create_data = folder_to_create_data
+
+ def rotation_matrix(self, input_vector, target_vector):
+ """
+ Computes the rotation matrix that aligns one 3D vector with another.
+
+ Args:
+ input_vector (array-like): The initial 3D vector.
+ target_vector (array-like): The target 3D vector.
+
+ Returns:
+ numpy.ndarray: A 3x3 rotation matrix.
+ """
+
+ # w_hat = np.cross(input_vector, target_vector)
+
+ theta = np.arccos(np.dot(input_vector, target_vector))
+
+ c = np.cos(theta)
+ s = np.sin(theta)
+ R = np.array([[1, 0, 0], [0, c, -s], [0, s, c]])
+
+ return R
+
+ def build(self, reference_position, aerodynamic_surface, include_control_devices=True):
+ """
+ Constructs the AVL-compatible aerodynamic surface from the given input data.
+
+ Args:
+ reference_position (dict): Reference position of the aerodynamic surface.
+ aerodynamic_surface (dict): Dictionary containing geometry and control device data.
+ include_control_devices (bool, optional): Whether to include control devices. Defaults to True.
+
+ Returns:
+ avl.Surface: The constructed AVL aerodynamic surface.
+ """
+ geometry = aerodynamic_surface["geometry"]
+ controls = aerodynamic_surface["control_devices"]
+
+ # Determine rotation matrix for geometry origin -> necessary due to extrusion of aerodynamic objects according to aircraftgeometry2
+ geometry_normal = np.array([geometry.normal.dx(), geometry.normal.dy(), geometry.normal.dz()])
+ extrusion_normal = np.array([0, 0, 1])
+ self.origin_unicado_to_avl = self.rotation_matrix(extrusion_normal, geometry_normal)
+
+ # Initialize avl sections list
+ avl_sections = []
+
+ # Set number of geometry sections
+ number_of_geometry_sections = len(geometry.sections)
+
+ for id in range(number_of_geometry_sections):
+
+ section = geometry.sections[id]
+ chord = section.get_chord_length()
+
+ # leading_edge_reference = np.dot(self.origin_unicado_to_avl, np.array(
+ # [section.origin.x(), section.origin.y(), section.origin.z()]))
+ leading_edge_reference = self.xyz_unicado_to_avl(
+ [section.origin.x(), section.origin.y(), section.origin.z()]) # leading_edge_reference.tolist()
+
+ # remove small elements from leading edge reference to avoid asymmetric behaviour
+ leading_edge_reference = [0.0 if abs(elem) < 1E-6 else elem for elem in leading_edge_reference]
+
+ scaling = section.get_thickness_scale()
+
+ # Profile scaling
+
+ # XFoil profile destination
+ xfoil_destination_dir = f"{self.folder_to_create_data}/profiles"
+ # Create xfoil_destintation_dir if not existing
+ os.makedirs(xfoil_destination_dir, exist_ok=True)
+ profile_path = f"{self.paths['airfoil_data_directory']}/{section.name}.dat"
+ xfoil_profile_path = f"{xfoil_destination_dir}/{section.name}.dat"
+
+ # Copy profile to xfoil_destination_dir -> needs to be copied due to limitation of profiles
+ # todo -> remove xfoil usage by scaling thickness directly
+ shutil.copy(profile_path, xfoil_profile_path)
+
+ section_profile_file_name = Profile(xfoil_profile_path, geometry.name,
+ "geometry", id).save_to_file(xfoil_destination_dir, scaling)
+
+ avl_section = avl.Section(refLE=leading_edge_reference, chord=chord,
+ dihedral=0.0, nspan=5, afileKey="AFILE", afile=section_profile_file_name)
+
+ avl_sections.append(avl_section)
+
+ if include_control_devices:
+ avl_sections = self.build_control_sections(avl_sections, geometry, controls)
+
+ self.surface = avl.Surface(tag=geometry.name)
+
+ for section in avl_sections:
+ self.surface.add_section(section)
+ # Do translation, angle etc at end
+ self.surface.add_translate(avl.Translate([reference_position["x"] + geometry.origin.x(),
+ reference_position["y"] + geometry.origin.y(),
+ reference_position["z"] + geometry.origin.z()]))
+ self.surface.add_angle(avl.Angle(geometry.rotation_z))
+ if geometry.is_symmetric: # symmetric
+ self.surface.add_yduplicate(avl.Yduplicate(0.0))
+
+ return self.surface
+
+ def xyz_unicado_to_avl(self, geom2_xyz=[]):
+ """
+ Converts coordinates from UNICADO format to AVL format.
+
+ Args:
+ geom2_xyz (list, optional): List of x, y, z coordinates in UNICADO format.
+
+ Returns:
+ list: Transformed coordinates in AVL format.
+ """
+ return np.dot(self.origin_unicado_to_avl, geom2_xyz).tolist()
+
+ def build_control_sections(self, avl_sections, geometry, controls):
+ """
+ Builds control surface sections for AVL.
+
+ Args:
+ avl_sections (list): List of existing AVL sections.
+ geometry (object): Aircraft geometry object.
+ controls (list): List of control surfaces.
+
+ Returns:
+ list: Updated AVL sections with control devices.
+ """
+ half_span = geom2.measure.span(geometry)
+ if geometry.is_symmetric:
+ half_span *= 0.5
+ sections_available = [section.origin.z() for section in geometry.sections]
+
+ # Set control surface vector parallel to surface normal
+ xyzhvec = [abs(geometry.normal.dx()), abs(geometry.normal.dy()), abs(geometry.normal.dz())]
+
+ for control in controls:
+ for idx, section in enumerate(control.sections):
+ # Build section
+ position_spanwise = abs(half_span * section.origin.z())*-1.0
+ current_id = None
+ for id, position in enumerate(sections_available):
+ if self.float_equality(position, position_spanwise):
+ current_id = id
+ break
+ else:
+ continue
+
+ # if current position is not existing -> create section
+ if current_id is None:
+ sections_available.append(position_spanwise)
+ origin = geom2.measure.offset_LE(geometry, position_spanwise)
+ leading_edge_reference = self.xyz_unicado_to_avl([origin.x(), origin.y(), origin.z()])
+ leading_edge_reference = [0.0 if abs(elem) < 1E-6 else elem for elem in leading_edge_reference]
+
+ id_left, id_right = self.find_adjacent_sections(geometry, position_spanwise)
+ position_scale, position_profile = self.find_thickness_scaling_and_profile(
+ geometry, id_left, id_right, position_spanwise)
+
+ chord = geom2.measure.chord(geometry, position_spanwise)
+
+ # Profile scaling
+ # XFoil profile destination
+ xfoil_destination_dir = f"{self.folder_to_create_data}/profiles"
+ # Create xfoil_destintation_dir if not existing
+ os.makedirs(xfoil_destination_dir, exist_ok=True)
+ profile_path = f"{self.paths['airfoil_data_directory']}/{position_profile}.dat"
+ xfoil_profile_path = f"{xfoil_destination_dir}/{section.name}.dat"
+
+ # Copy profile to xfoil_destination_dir -> needs to be copied due to limitation of profiles
+ # todo -> remove xfoil usage by scaling thickness directly
+ shutil.copy(profile_path, xfoil_profile_path)
+
+ section_profile_file_name = Profile(xfoil_profile_path, geometry.name,
+ control.name, idx).save_to_file(xfoil_destination_dir, position_scale)
+
+ avl_section = avl.Section(refLE=leading_edge_reference, chord=chord,
+ dihedral=0.0, nspan=5, afileKey="AFILE", afile=section_profile_file_name)
+
+ contour = section.get_contour(False)
+
+
+ # Hinge line according to avl documentation
+ # LE Surface 0 ... -xhinge
+ # TE Surface xhinge ... 1
+
+ # Assume trailing edge device
+ xhinge = contour.vertex(0).x()
+ # Change if xhinge is at zero
+ if abs(xhinge) < 1E-4:
+ xhinge = -contour.vertex(1).x()
+ # Direction for control surface - according to avl documentation
+ sgndup = 1
+ gain = 1
+ if control.name.startswith("aileron"):
+ sgndup = -1
+ if control.name.startswith("rudder") or control.name.startswith("slat"):
+ gain = -1
+
+ # Create avl control device
+ control_device = avl.Control(tag=control.name, gain=gain, xyzhvec=xyzhvec,
+ xhinge=xhinge, sgndup=sgndup)
+
+ # Add control device section
+ avl_section.add_control_device(control_device)
+ avl_sections.append(avl_section)
+ else:
+ # add_control to section
+ contour = section.get_contour(False)
+ # Hinge line according to avl documentation
+ # LE Surface 0 ... -xhinge
+ # TE Surface xhinge ... 1
+
+ # Assume trailing edge device
+ xhinge = contour.vertex(0).x()
+ # Change if xhinge is at zero
+ if abs(xhinge) < 1E-4:
+ xhinge = -contour.vertex(1).x()
+ # Direction for control surface - according to avl documentation
+ sgndup = 1
+ gain = 1
+ if control.name.startswith("aileron_") or control.name == "aileron":
+ sgndup = -1
+ if control.name.startswith("rudder") or control.name.startswith("slat"):
+ gain = -1
+
+ control_device = avl.Control(tag=control.name, gain=gain, xyzhvec=xyzhvec,
+ xhinge=xhinge, sgndup=sgndup)
+ avl_sections[current_id].add_control_device(control_device)
+
+ # Sort avl_sections
+ sorted_section_idx = sorted(range(len(sections_available)), key=lambda x: sections_available[x], reverse=True)
+ reordered_sections = [avl_sections[i] for i in sorted_section_idx]
+
+ # Get control section ids start + end
+ control_devices = defaultdict(list)
+ for idx, section in enumerate(reordered_sections):
+ for control in section.controls:
+ control_devices[control.tag].append(idx)
+
+ # Add control to section if control sections not next to eachother
+ for device, boundary_ids in control_devices.items():
+ first_section_id = boundary_ids[0]
+ last_section_id = boundary_ids[-1]
+ missing_section_ids = list(range(first_section_id+1,last_section_id))
+
+ # Append missing control to existing section
+ boundary_controls = [control for controls in [reordered_sections[first_section_id].controls, reordered_sections[last_section_id].controls] for control in controls if control.tag == device]
+ for id in missing_section_ids:
+ # Get first and last control devices
+ first_hinge = boundary_controls[0].xhinge
+ first_gain = boundary_controls[0].gain
+ last_hinge = boundary_controls[-1].xhinge
+ last_gain = boundary_controls[-1].gain
+
+ control_to_add = boundary_controls[0]
+ if not self.float_equality(first_hinge, last_hinge):
+ # Get first and last position in main extrusion direction -> xyzhvec
+ first_pos = [reordered_sections[first_section_id].ref[i] for i, val in enumerate(xyzhvec) if abs(val) > 0.9][0]
+ last_pos = [reordered_sections[last_section_id].ref[i] for i, val in enumerate(xyzhvec) if abs(val) > 0.9][0]
+
+ current_pos = [reordered_sections[id].ref[i] for i, val in enumerate(xyzhvec) if abs(val) > 0.9][0]
+
+ t = (current_pos-first_pos)/(last_pos-first_pos)
+ current_xhinge = self.lerp(first_hinge, last_hinge, t)
+ current_gain = self.lerp(first_gain,last_gain, t)
+ control_to_add.xhinge = current_xhinge
+ control_to_add.gain = current_gain
+
+ reordered_sections[id].controls.append(control_to_add)
+
+ # Assume that a device is always unique (e.g. slats are not separated -> otherwise it would be slat_1, slat_2 etc. )
+
+
+
+ return reordered_sections
+
+ def lerp(self, a: float, b: float, t: float) -> float:
+ """
+ Performs linear interpolation between two values.
+
+ Args:
+ a (float): The starting value.
+ b (float): The ending value.
+ t (float): The interpolation factor (0 ≤ t ≤ 1).
+
+ Returns:
+ float: The interpolated value.
+ """
+ return (1-t)* a + t*b
+
+ def float_equality(self, value, to_check, tol=1e-4):
+ """
+ Checks if two floating-point numbers are equal within a given tolerance.
+
+ Args:
+ value (float): The first number to compare.
+ to_check (float): The second number to compare.
+ tol (float, optional): The allowed tolerance for equality. Defaults to 1e-4.
+
+ Returns:
+ bool: True if values are considered equal, otherwise False.
+ """
+ return abs(value-to_check) < tol
+
+
+ def find_adjacent_sections(self, geometry, position):
+ """
+ Finds the indices of the two adjacent sections in the given geometry closest to the specified position.
+
+ Args:
+ geometry (object): The geometry containing sections.
+ position (float): The spanwise position to find the nearest sections.
+
+ Returns:
+ tuple[int, int]: Indices of the left and right adjacent sections.
+ """
+ section_positions = [section.origin.z() for section in geometry.sections]
+ id_left = 0
+ id_right = 1
+ if not position in section_positions:
+ for idx, section_position in enumerate(section_positions):
+ if position <= section_position:
+ id_left = idx
+ id_right = idx+1
+
+ return id_left, id_right
+
+ def find_thickness_scaling_and_profile(self, geometry, left_section_id, right_section_id, position):
+ """
+ Determines the thickness scaling and profile name based on the given spanwise position.
+
+ Args:
+ geometry (object): The aerodynamic geometry.
+ left_section_id (int): The index of the left section.
+ right_section_id (int): The index of the right section.
+ position (float): The spanwise position to evaluate.
+
+ Returns:
+ tuple[float, str]: The interpolated thickness scaling factor and the selected airfoil profile name.
+ """
+ left_section = geometry.sections[left_section_id]
+ left_section_profile = left_section.name
+ left_section_scale = left_section.get_thickness_scale()
+
+ right_section = geometry.sections[right_section_id]
+ right_section_profile = right_section.name
+ right_section_scale = right_section.get_thickness_scale()
+
+ normalized_position = (position-left_section.origin.z())/(right_section.origin.z()-left_section.origin.z())
+
+ position_scale: float = left_section_scale * \
+ (1 - normalized_position) + right_section_scale * normalized_position
+ position_profile: str = left_section_profile
+ if right_section_profile != left_section_profile:
+ if normalized_position > 0.5:
+ position_profile = right_section_profile
+ else:
+ position_profile = left_section_profile
+
+ return position_scale, position_profile
+
+ # left_section_profile
+
+
+
+class Profile:
+ """Profile storage and scaling class."""
+ def __init__(self, profile_file: str, section_geometry_name: str, section_type: str = "geometry", section_id: int = 0) -> None:
+ """
+ Initializes the Profile class.
+
+ Args:
+ profile_file (str): The filename of the profile.
+ section_geometry_name (str): The name of the geometry section.
+ section_type (str, optional): The type of the section. Defaults to "geometry".
+ section_id (int, optional): The ID of the section. Defaults to 0.
+ """
+
+ self.profile_file = profile_file
+ self.section_id = section_id
+ self.output_file = f"section_{section_geometry_name}_{section_type}_{section_id}.dat"
+
+ def remove_existing_profile(self, path: str) -> None:
+ """
+ Removes an existing profile file if it exists.
+
+ Args:
+ path (str): The directory path where the profile file is located.
+ """
+ path_to_file = f"{path}/{self.output_file}"
+ if os.path.exists(path_to_file):
+ os.remove(path_to_file)
+ else:
+ pass
+
+ def save_to_file(self, path: str, scale: float = 1.) -> str:
+ """
+ Applies thickness scaling to the profile and saves the modified profile.
+
+ Args:
+ path (str, optional): The directory path for saving the scaled profile. Defaults to an empty string.
+ scale (float, optional): The scaling factor for thickness. Defaults to 1.0.
+ """
+
+ self.remove_existing_profile(path)
+
+ # Scale thickness
+ self.run_thickness_scaling(path, scale)
+ return f"{path}/{self.output_file}"
+
+ def run_thickness_scaling(self, path:str="", scale: float=1.0) -> None:
+ """ Runs thickness scaling
+
+ Args:
+ path (str, optional): path to airfoil directory. Defaults to "".
+ scale (float, optional): scales thickness according to value. Defaults to 1.0.
+ """
+ output_path = f"{path}/{self.output_file}"
+
+ coordinates = load_profile(self.profile_file)
+ scaled_coordinates = rescale_profile(coordinates=coordinates, scale_t_to_c=scale)
+
+ save_profile(scaled_coordinates, output_path)
diff --git a/pyavlunicadopackage/pyavlunicado/avlunicadoio.py b/pyavlunicadopackage/pyavlunicado/avlunicadoio.py
new file mode 100644
index 0000000000000000000000000000000000000000..f0e808ecd6e4fc9fde507efcd4c8a4494d328a1d
--- /dev/null
+++ b/pyavlunicadopackage/pyavlunicado/avlunicadoio.py
@@ -0,0 +1,290 @@
+# UNICADO - UNIversity Conceptual Aircraft Design and Optimization
+#
+# Copyright (C) 2025 UNICADO consortium
+#
+# 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/>.
+#
+# Description:
+# This file is part of UNICADO.
+
+"""
+PyAvlUnicado - A Python Interface for Athena Vortex Lattice (AVL) for use with UNICADO
+
+This module provides a Python interface for generating and executing AVL simulations. It includes:
+- Geometry initialization and manipulation
+- Run case setup and execution
+- Interaction with AVL for aerodynamic analysis
+
+Features:
+- Reads and processes aircraft geometry
+- Manages aerodynamic surfaces and control devices
+- Creates AVL input files for simulations
+- Automates AVL run case execution and result extraction
+
+Dependencies:
+- numpy
+- ambiance
+- pyaircraftgeometry2
+- pyavl
+"""
+import pyaircraftgeometry2 as geom2
+import pyavl as avl
+import os
+import shutil
+from .avlunicadogeometry import AerodynamicSurface
+from .avlunicadoutility import load_aerodynamic_surfaces
+from ambiance import Atmosphere
+import numpy as np
+
+
+class Avlinterface:
+ """
+ Manages the interface for AVL simulations, including geometry initialization and run case management.
+ """
+ aerodynamic_surfaces = None
+ aerodynamic_components: list=[]
+ aerodynamic_reference_data: list=[]
+ aerodynamic_surface_control_keys = set()
+ paths: str=None
+ aircraft_name: str=None
+ runcases: list=[]
+ derivative_data = None
+ none_high_lift_devices: list = ["aileron", "elevator", "rudder"]
+
+ def __init__(self, paths: str = None, aircraft_name: str = "default_ac") -> None:
+ """
+ Initializes the Avlinterface class.
+
+ Args:
+ paths (str, optional): Path to aircraft geometry files.
+ aircraft_name (str, optional): Name of the aircraft model.
+ """
+ self.paths = paths
+ self.aircraft_name = aircraft_name
+
+ def initialize_geometry(self, component_list: list[str] = [], use_control_devices=False, use_high_lift_devices=False):
+ """
+ Initializes aerodynamic geometry based on input components.
+
+ Args:
+ component_list (list, optional): List of components to include.
+ use_control_devices (bool, optional): Whether to include control devices.
+ use_high_lift_devices (bool, optional): Whether to include high lift devices.
+ """
+ # clean up aerodynamic components
+ if self.aerodynamic_components:
+ self.aerodynamic_components.clear()
+
+ # Load aerodynamic surfaces from acxml
+ self.aerodynamic_surfaces = load_aerodynamic_surfaces(self.paths, component_list)
+
+ # Read reference data (S_ref, mac, Span)
+ if "wing" in self.aerodynamic_surfaces:
+ reference_area: float = 0
+ reference_mac: float = 0
+ reference_span: float = 0
+ for aerodynamic_surface in self.aerodynamic_surfaces["wing"]["aerodynamic_surfaces"]:
+ current_wing = aerodynamic_surface["geometry"]
+ current_area = geom2.measure.reference_area(current_wing)
+
+ # select reference values from biggest area
+ if current_area > reference_area:
+ reference_area = current_area
+ reference_mac = geom2.measure.mean_aerodynamic_chord(current_wing)
+ reference_span = geom2.measure.span(current_wing)
+ reference_phi_25 = -geom2.measure.sweep(current_wing, 0.25*reference_span, 0.25)
+ self.aerodynamic_reference_data = [reference_area, reference_mac, reference_span, reference_phi_25]
+
+ for component in self.aerodynamic_surfaces:
+ # Read reference position
+ component_reference_position = self.aerodynamic_surfaces[component]["reference_position"]
+ for aerodynamic_surface in self.aerodynamic_surfaces[component]["aerodynamic_surfaces"]:
+
+ # If use_high_lift_devices -> False ... only none_high_lift_devices are used for geometry, else all elements are used
+ # WARNING -> if all devices are used, make sure that they are not too close to eachother, otherwise avl has issues with computation
+ if not use_high_lift_devices:
+ aerodynamic_surface["control_devices"] = [control for control in aerodynamic_surface["control_devices"] if control.name in self.none_high_lift_devices]
+
+ self.aerodynamic_components.append(AerodynamicSurface(self.paths, self.aircraft_name).build(
+ component_reference_position, aerodynamic_surface, use_control_devices))
+
+
+ for control in aerodynamic_surface["control_devices"]:
+ self.aerodynamic_surface_control_keys.add(control.name)
+
+
+ def create_avl_geometry_file(self, author: str = "", filename: str = "", runcasename: str = "default", mach_number: float = 0.78, reference_data: list[float] = None, reference_cog: list[float] = None):
+ """
+ Creates an AVL geometry input file.
+ """
+ if filename == "":
+ filename = f"{self.aircraft_name}/{self.aircraft_name}"
+ else:
+ filename = f"{self.aircraft_name}/{filename}"
+ if reference_data is None:
+ reference_data = self.aerodynamic_reference_data
+ if reference_cog is None:
+ reference_cog = [0., 0., 0.]
+ print(f"Reference center of gravity is None ... set to {reference_cog}")
+
+ header_avl_file = [avl.Header(runcase=runcasename, mach=mach_number,
+ sym=[0, 0, 0.0], base=reference_data[:3], ref=reference_cog)]
+
+ avl.write_avl_file(author, header_avl_file + self.aerodynamic_components, filename)
+
+
+
+ def add_runcase(self, runcase):
+ """ Add runcase
+
+ Args:
+ runcase (object): Runcase
+ """
+ if not isinstance(runcase, AvlRuncase):
+ print("Runcase will not be appended - no AvlRuncase...")
+
+ else:
+ self.runcases.append(runcase)
+
+ def number_of_runcases(self):
+ """Return numer of runcases"""
+ return len(self.runcases)
+
+ def runcase_ids(self):
+ """Return list of runcase ids"""
+ return [rc.runcase_id for rc in self.runcases]
+
+ def define_runcase_batch(self, runcase_folder: str=".", alphas: list[float]=[0.0], betas: list[float]=[0.0], pb2vs: list[float]=[0.0], qc2vs: list[float]=[0.0], rb2vs: list[float]=[0.0], mach_numbers: list[float]=[0.0], CD0: float=0.02, heights_in_m: list[float]=[0.0], cog: list[float]=[0.0,0.0,0.0], controls: list[str]=[], clear_runcase_folder: bool=True):
+ """
+ Defines a batch of AVL run cases by varying aerodynamic parameters.
+
+ Args:
+ runcase_folder (str, optional): Path to store run case files. Defaults to ".".
+ alphas (list[float], optional): List of angle of attack values. Defaults to [0.0].
+ betas (list[float], optional): List of sideslip angle values. Defaults to [0.0].
+ pb2vs (list[float], optional): List of roll rate values. Defaults to [0.0].
+ qc2vs (list[float], optional): List of pitch rate values. Defaults to [0.0].
+ rb2vs (list[float], optional): List of yaw rate values. Defaults to [0.0].
+ mach_numbers (list[float], optional): List of Mach numbers. Defaults to [0.0].
+ CD0 (float, optional): Zero-lift drag coefficient. Defaults to 0.02.
+ heights_in_m (list[float], optional): List of altitude values in meters. Defaults to [0.0].
+ cog (list[float], optional): Center of gravity coordinates [x, y, z]. Defaults to [0.0, 0.0, 0.0].
+ controls (list[str], optional): List of control surface names. Defaults to an empty list.
+ clear_runcase_folder (bool, optional): Whether to clear the run case folder before creation. Defaults to True.
+ """
+ runcase_id = self.number_of_runcases() + 1
+ print("Clear existing runcases ...")
+ self.runcases.clear()
+ if os.path.isdir(runcase_folder) and clear_runcase_folder:
+ if not self.is_directory_empty(runcase_folder):
+ self.empty_directory(runcase_folder)
+
+ # Create list from set
+ if len(controls) == 0:
+ controls=list(self.aerodynamic_surface_control_keys)
+
+ # Create different runcase batch
+ for mach_number, height_in_m in zip(mach_numbers,heights_in_m):
+ for rb2v in rb2vs:
+ for qc2v in qc2vs:
+ for pb2v in pb2vs:
+ for beta in betas:
+ for alpha in alphas:
+ self.runcases.append(AvlRuncase(runcase_folder, runcase_id, f"RunCase {runcase_id} - {height_in_m:.2f}", alpha, beta, pb2v, qc2v, rb2v, mach_number*np.cos(self.aerodynamic_reference_data[3]), CD0, height_in_m, cog, controls).create())
+ runcase_id += 1
+
+ print(f"Current number of Runcases: ... {self.number_of_runcases()}")
+
+ def run_avl(self, runcase_folder, stability_files_folder, results_folder, derivative_collection_file=None):
+ """
+ Executes an AVL simulation using the specified run cases and stability file storage.
+
+ Args:
+ runcase_folder (str): Path to the folder containing AVL run case files.
+ stability_files_folder (str): Path to the folder where AVL stability results will be stored.
+ results_folder (str): Path to the folder where final results will be saved.
+ derivative_collection_file (str, optional): Name of the output CSV file to store aerodynamic derivatives.
+
+ Returns:
+ avl.PyAvlStabilityFilesConvert or None: Returns processed aerodynamic derivatives if `derivative_collection_file` is provided, otherwise None.
+ """
+ geometry_file_name = f"{self.aircraft_name}.avl"
+ runner = avl.pyavlrunner.PyAvlRunner(geometry_file_name)
+ runner.session_setup_single_rc(path_to_session_folder=self.aircraft_name, runcase_files_folder=runcase_folder, stability_files_folder=stability_files_folder, results_folder=results_folder, runcase_ids=self.runcase_ids())
+ runner.run_avl_session()
+
+ if isinstance(derivative_collection_file, str):
+
+ os.makedirs(results_folder, exist_ok=True)
+ stability_files = [f"{runner.stability_files_folder()}/{file}" for file in os.listdir(runner.stability_files_folder()) if file.endswith(".stab")]
+ self.derivative_data = avl.PyAvlStabilityFilesConvert(filenames=stability_files,
+ keyListControls=list(self.aerodynamic_surface_control_keys), fileout=f"{results_folder}/{derivative_collection_file}.csv")
+ return self.derivative_data
+
+ def is_directory_empty(self,directory):
+ """Check whether a directory is empty or ot"""
+ return not any(os.scandir(directory))
+
+ def empty_directory(self,directory):
+ """ Empty directory
+
+ Args:
+ directory (str): removes the content of this directory
+ """
+ for item in os.listdir(directory):
+ item_path = os.path.join(directory, item)
+ if os.path.isfile(item_path):
+ os.remove(item_path)
+ elif os.path.isdir(item_path):
+ shutil.rmtree(item_path)
+
+class AvlRuncase:
+ """
+ Defines an AVL run case with aerodynamic parameters.
+ """
+ def __init__(self, runcase_folder: str=".", runcase_id: int=1, runcase_name: str="default_runcase", alpha: float=0.0, beta: float=0.0, pb2v: float=0.0, qc2v: float=0.0, rb2v: float=0.0, mach_number: float=0.0, CD0: float=0.02, height_in_m: float=0.0, cog: list[float]=[0.0,0.0,0.0], controls=None):
+ """
+ Initializes an AVL run case.
+ """
+ self.runcase_directory = runcase_folder
+ self.runcase_id = runcase_id
+ self.runcase_name = runcase_name
+ self.runcase_file: str= f"rc_{runcase_id:03}.run"
+ self.alpha = alpha
+ self.beta = beta
+ self.pb2v = pb2v
+ self.qc2v = qc2v
+ self.rb2v = rb2v
+ self.controls = controls
+ self.CD0 = CD0
+ self.mach_number = mach_number
+ self.altitude = height_in_m
+ self.density = Atmosphere(self.altitude).density
+ self.speed_of_sound = Atmosphere(self.altitude).speed_of_sound
+ self.velocity = mach_number * self.speed_of_sound
+ self.cog = cog
+
+ def create(self):
+ """
+ Creates an AVL run case file in the specified directory.
+ """
+ # Create runcase directory
+ os.makedirs(self.runcase_directory, exist_ok=True)
+
+ full_path_to_runcase = f"{self.runcase_directory}/{self.runcase_file}"
+ if os.path.exists(full_path_to_runcase):
+ os.remove(full_path_to_runcase)
+ with open(full_path_to_runcase, "+w") as rcfile:
+ avl.write_runcase(file=rcfile, runcasenum=1, runcasename=self.runcase_name, alpha=self.alpha, beta=self.beta, pb2v=self.pb2v, qc2v=self.qc2v, rb2v=self.rb2v, mach=self.mach_number, CD0=self.CD0, density=self.density, a=self.speed_of_sound, cgref=self.cog, controls=self.controls)
+ return self
diff --git a/pyavlunicadopackage/pyavlunicado/avlunicadoutility.py b/pyavlunicadopackage/pyavlunicado/avlunicadoutility.py
new file mode 100644
index 0000000000000000000000000000000000000000..bc7a79277726c3e9ebe57904946228aeb8a55bc2
--- /dev/null
+++ b/pyavlunicadopackage/pyavlunicado/avlunicadoutility.py
@@ -0,0 +1,179 @@
+# UNICADO - UNIversity Conceptual Aircraft Design and Optimization
+#
+# Copyright (C) 2025 UNICADO consortium
+#
+# 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/>.
+#
+# Description:
+# This file is part of UNICADO.
+
+"""
+PyAvlUnicado - A Python Interface for Athena Vortex Lattice (AVL) for use with UNICADO
+
+This module provides utility functions for handling aircraft geometry data,
+specifically for interacting with aircraft XML files and extracting aerodynamic
+surface information.
+
+Functions:
+ - get_element_path(element, tree): Retrieves the hierarchical path of an XML element.
+ - load_aerodynamic_surfaces(paths, component_names): Loads aerodynamic surfaces from aircraft XML data.
+ - load_aerodynamic_surfaces_from_information(paths, aerodynamic_surfaces_information): Extracts aerodynamic surface details.
+ - aerodynamic_surfaces_from_component_available(component): Identifies aerodynamic surfaces in a component.
+ - load_reference_position(component): Extracts the reference position of a component.
+"""
+
+import pyaircraftgeometry2 as geom2
+import pyaixml as aixml
+
+
+def get_element_path(element, tree):
+ """
+ Retrieves the hierarchical path of an XML element relative to the root.
+
+ Args:
+ element (ET.Element): The XML element whose path is to be determined.
+ tree (ET.Element): The root of the XML tree.
+
+ Returns:
+ str: The hierarchical path of the element.
+ """
+ element_path = []
+ while element != tree:
+ parent = tree.find(f".//{element.tag}/..")
+ element_path.append(element.tag)
+ element = parent
+ element_path.reverse()
+ return '/'.join(element_path)
+
+
+def load_aerodynamic_surfaces(paths, component_names: list[str]=[]):
+ """
+ Loads aerodynamic surface data from aircraft XML files.
+
+ Args:
+ paths (dict): Dictionary containing paths to aircraft exchange files.
+ component_names (list[str], optional): List of component names to extract data from.
+
+ Returns:
+ dict: A dictionary containing aerodynamic surface details for each component.
+
+ Raises:
+ RuntimeError: If component_names is not a valid list of strings.
+ """
+ component_design_items = None
+ if component_names or isinstance(component_names, list) and all(isinstance(name, str) for name in component_names):
+
+ if component_names:
+ component_design_items = []
+ for component_name in component_names:
+ component_design_items.append(paths["root_of_aircraft_exchange_tree"].find(
+ f".//component_design/{component_name}"))
+
+ else:
+ component_design_items = paths["root_of_aircraft_exchange_tree"].find('.//component_design')
+
+ aerodynamic_surfaces = {}
+ if component_design_items is not None:
+ for component in component_design_items:
+ info = aerodynamic_surfaces_from_component_available(component)
+ if info is not None:
+ aerodynamic_surfaces[component.tag] = {
+ "aerodynamic_surfaces": load_aerodynamic_surfaces_from_information(paths, info),
+ "reference_position": load_reference_position(component),
+ "info": info
+ }
+ return aerodynamic_surfaces
+ else:
+ raise RuntimeError(
+ f"component_names does not match required type list - current type: {type(component_names)}")
+
+
+def load_aerodynamic_surfaces_from_information(paths, aerodynamic_surfaces_information):
+ """
+ Extracts aerodynamic surfaces from the given aircraft XML data.
+
+ Args:
+ paths (dict): Dictionary containing paths to aircraft exchange files.
+ aerodynamic_surfaces_information (list): List of aerodynamic surface information.
+
+ Returns:
+ list: A list of aerodynamic surfaces including geometry and control devices.
+ """
+ acxml = aixml.openDocument(paths["path_to_aircraft_exchange_file"])
+ aerodynamic_surfaces = []
+ for aerodynamic_surface_info in aerodynamic_surfaces_information:
+ aerodynamic_surface = {}
+ aerodynamic_surface["geometry"] = geom2.factory.WingFactory(
+ acxml, paths["airfoil_data_directory"]).create(aerodynamic_surface_info["geometry_path"])
+
+ aerodynamic_surface["control_devices"] = []
+ for control_device in aerodynamic_surface_info["control_devices"]:
+ device = geom2.factory.ControlDeviceFactory(acxml, paths["airfoil_data_directory"]).create(control_device["path"])
+ device.name = control_device["name"]
+ aerodynamic_surface["control_devices"].append(device)
+ aerodynamic_surfaces.append(aerodynamic_surface)
+ return aerodynamic_surfaces
+
+
+def aerodynamic_surfaces_from_component_available(component):
+ """
+ Identifies aerodynamic surfaces in an aircraft component.
+
+ Args:
+ component (ET.Element): XML element representing the aircraft component.
+
+ Returns:
+ list or None: A list of aerodynamic surfaces found in the component, or None if no surfaces are available.
+ """
+ components_aerodynamic_surfaces = component.findall(".//aerodynamic_surface[@ID]")
+ if components_aerodynamic_surfaces:
+ aerodynamic_surfaces_available = []
+
+ for components_aerodynamic_surface in components_aerodynamic_surfaces:
+ aerodynamic_surface_available = {}
+ aerodynamic_surface_available["geometry_path"] = f"{component.tag}/{get_element_path(components_aerodynamic_surface, component)}@{components_aerodynamic_surface.attrib['ID']}"
+ aerodynamic_surface_available["control_devices"] = []
+
+ # look for control devices
+ control_devices = components_aerodynamic_surface.findall(".//control_device[@ID]")
+ if control_devices:
+ for control_device in control_devices:
+ control_device_info = {
+ "path": f"{aerodynamic_surface_available['geometry_path']}/{get_element_path(control_device, components_aerodynamic_surface)}@{control_device.attrib['ID']}",
+ "name": control_device.attrib["description"]
+ }
+ aerodynamic_surface_available["control_devices"].append(control_device_info)
+ aerodynamic_surfaces_available.append(aerodynamic_surface_available)
+ return aerodynamic_surfaces_available
+ else:
+ return None
+
+def load_reference_position(component):
+ """
+ Extracts the reference position of an aircraft component.
+
+ Args:
+ component (ET.Element): XML element representing the aircraft component.
+
+ Returns:
+ dict: Dictionary containing x, y, and z coordinates of the reference position.
+ """
+ position = {
+ "x": float(component.find("./position/x/value").text),
+ "y": float(component.find("./position/y/value").text),
+ "z": float(component.find("./position/z/value").text)
+ }
+
+ return position
+
diff --git a/pyavlunicadopackage/pyavlunicado/xfoilgeometryscale.py b/pyavlunicadopackage/pyavlunicado/xfoilgeometryscale.py
new file mode 100644
index 0000000000000000000000000000000000000000..0f9d1df2ba69c5391ae902e521e964c342782b9b
--- /dev/null
+++ b/pyavlunicadopackage/pyavlunicado/xfoilgeometryscale.py
@@ -0,0 +1,885 @@
+# UNICADO - UNIversity Conceptual Aircraft Design and Optimization
+#
+# Copyright (C) 2025 UNICADO consortium
+#
+# 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/>.
+#
+# Description:
+# This file is part of UNICADO.
+
+"""
+------------------------------------------------------------------------------
+This file is a direct translation and copy of portions of the XFOIL code,
+originally written by Mark Drela. The implementation below is based on
+the XFOIL routines and is provided as-is under the original XFOIL license.
+
+XFOIL License:
+--------------
+Copyright (C) 1996, 1998, 2002, 2005 by Mark Drela. All rights reserved.
+
+XFOIL is provided "as is" without warranty of any kind, either expressed
+or implied, including but not limited to the implied warranties of
+merchantability and fitness for a particular purpose.
+
+For more information on the XFOIL license, please refer to the license file
+included with the original XFOIL distribution or visit:
+http://web.mit.edu/drela/Public/web/xfoil/
+
+------------------------------------------------------------------------------
+"""
+
+"""
+xfoil_geometry.py
+
+A collection of routines to compute spline derivatives,
+evaluate splines, find leading edge parameters, opposite points,
+arc lengths, curvatures, geometric properties, and to update airfoil
+camber and thickness.
+"""
+
+import numpy as np
+
+
+def tri_solve(A: np.ndarray, B: np.ndarray, C: np.ndarray, D: np.ndarray, kk: int) -> None:
+ """
+ Solve a tri-diagonal system in place.
+
+ The system is assumed to be:
+ A[i]*D[i] + B[i+1]*D[i+1] + C[i]*D[i-1] = original_D[i]
+ where D is replaced with the solution.
+
+ Parameters
+ ----------
+ A : ndarray
+ Main diagonal (length kk).
+ B : ndarray
+ Upper diagonal (length kk); note B[0] is unused.
+ C : ndarray
+ Lower diagonal (length kk); note C[-1] is unused.
+ D : ndarray
+ Right-hand side (length kk); will be overwritten with solution.
+ kk : int
+ Number of equations.
+ """
+ for k in range(1, kk):
+ km = k - 1
+ C[km] /= A[km]
+ D[km] /= A[km]
+ A[k] -= B[k] * C[km]
+ D[k] -= B[k] * D[km]
+ D[kk - 1] /= A[kk - 1]
+ for k in range(kk - 2, -1, -1):
+ D[k] -= C[k] * D[k + 1]
+ # The solution is now in D.
+
+
+def splind(X: np.ndarray, S: np.ndarray, xs1: float, xs2: float) -> np.ndarray:
+ """
+ Compute spline derivative array for a function defined by X(S).
+
+ End conditions:
+ - If xs1 or xs2 equal 999.0, use the usual zero second derivative condition.
+ - If xs1 or xs2 equal -999.0, use the zero third derivative condition.
+ - Otherwise, use the provided derivative.
+
+ Parameters
+ ----------
+ X : ndarray
+ Dependent variable array.
+ S : ndarray
+ Independent variable (parameter) array.
+ xs1 : float
+ Left endpoint derivative condition.
+ xs2 : float
+ Right endpoint derivative condition.
+
+ Returns
+ -------
+ ndarray
+ Spline derivative array.
+ """
+ N = len(X)
+ if N > 1000:
+ raise ValueError("SPLIND: array overflow, increase NMAX")
+ A = np.zeros(N)
+ B = np.zeros(N)
+ C = np.zeros(N)
+ D = np.zeros(N)
+
+ # Interior points (i = 1 to N-2)
+ for i in range(1, N - 1):
+ dsm = S[i] - S[i - 1]
+ dsp = S[i + 1] - S[i]
+ B[i] = dsp
+ A[i] = 2.0 * (dsm + dsp)
+ C[i] = dsm
+ D[i] = 3.0 * (((X[i + 1] - X[i]) * dsm / dsp) +
+ ((X[i] - X[i - 1]) * dsp / dsm))
+ # Left endpoint condition
+ if xs1 == 999.0:
+ A[0] = 2.0
+ C[0] = 1.0
+ D[0] = 3.0 * (X[1] - X[0]) / (S[1] - S[0])
+ elif xs1 == -999.0:
+ A[0] = 1.0
+ C[0] = 1.0
+ D[0] = 2.0 * (X[1] - X[0]) / (S[1] - S[0])
+ else:
+ A[0] = 1.0
+ C[0] = 0.0
+ D[0] = xs1
+
+ # Right endpoint condition
+ if xs2 == 999.0:
+ B[-1] = 1.0
+ A[-1] = 2.0
+ D[-1] = 3.0 * (X[-1] - X[-2]) / (S[-1] - S[-2])
+ elif xs2 == -999.0:
+ B[-1] = 1.0
+ A[-1] = 1.0
+ D[-1] = 2.0 * (X[-1] - X[-2]) / (S[-1] - S[-2])
+ else:
+ A[-1] = 1.0
+ B[-1] = 0.0
+ D[-1] = xs2
+ if N == 2 and xs1 == -999.0 and xs2 == -999.0:
+ B[-1] = 1.0
+ A[-1] = 2.0
+ D[-1] = 3.0 * (X[1] - X[0]) / (S[1] - S[0])
+ tri_solve(A, B, C, D, N)
+ return D
+
+
+def segspl(X: np.ndarray, S: np.ndarray) -> np.ndarray:
+ """
+ Compute segmented spline derivatives for X(S), allowing discontinuities at duplicate S values.
+
+ Parameters
+ ----------
+ X : ndarray
+ Data array.
+ S : ndarray
+ Parameter array.
+
+ Returns
+ -------
+ ndarray
+ Derivative array computed piecewise using splind.
+
+ Raises
+ ------
+ ValueError
+ If the first or last S values are duplicated.
+ """
+ N = len(S)
+ if np.isclose(S[0], S[1]):
+ raise ValueError("SEGSPL: First input point duplicated")
+ if np.isclose(S[-1], S[-2]):
+ raise ValueError("SEGSPL: Last input point duplicated")
+ XS = np.empty_like(X)
+ seg_start = 0
+ for i in range(1, N - 1):
+ if np.isclose(S[i], S[i + 1]):
+ seg_slice = slice(seg_start, i + 1)
+ XS[seg_slice] = splind(X[seg_slice].copy(), S[seg_slice].copy(), xs1=-999.0, xs2=-999.0)
+ seg_start = i + 1
+ seg_slice = slice(seg_start, N)
+ XS[seg_slice] = splind(X[seg_slice].copy(), S[seg_slice].copy(), xs1=-999.0, xs2=-999.0)
+ return XS
+
+
+def seval(ss: float, X: np.ndarray, XS: np.ndarray, S: np.ndarray) -> float:
+ """
+ Evaluate the spline defined by X(S) with derivatives XS at parameter value ss.
+
+ Parameters
+ ----------
+ ss : float
+ Parameter value at which to evaluate.
+ X : ndarray
+ Data array.
+ XS : ndarray
+ Spline derivative array.
+ S : ndarray
+ Parameter array.
+
+ Returns
+ -------
+ float
+ Interpolated value.
+ """
+ i = np.searchsorted(S, ss, side='right')
+ if i == 0:
+ i = 1
+ if i >= len(S):
+ i = len(S) - 1
+ ds = S[i] - S[i - 1]
+ t = (ss - S[i - 1]) / ds
+ cx1 = ds * XS[i - 1] - X[i] + X[i - 1]
+ cx2 = ds * XS[i] - X[i] + X[i - 1]
+ return t * X[i] + (1.0 - t) * X[i - 1] + (t - t * t) * ((1.0 - t) * cx1 - t * cx2)
+
+
+def deval(ss: float, X: np.ndarray, XS: np.ndarray, S: np.ndarray) -> float:
+ """
+ Evaluate the first derivative of the spline at parameter value ss.
+
+ Parameters
+ ----------
+ ss : float
+ Parameter value.
+ X : ndarray
+ Data array.
+ XS : ndarray
+ Spline derivative array.
+ S : ndarray
+ Parameter array.
+
+ Returns
+ -------
+ float
+ First derivative dX/dS evaluated at ss.
+ """
+ i = np.searchsorted(S, ss, side='right')
+ if i == 0:
+ i = 1
+ if i >= len(S):
+ i = len(S) - 1
+ ds = S[i] - S[i - 1]
+ t = (ss - S[i - 1]) / ds
+ cx1 = ds * XS[i - 1] - X[i] + X[i - 1]
+ cx2 = ds * XS[i] - X[i] + X[i - 1]
+ return (X[i] - X[i - 1] +
+ (1.0 - 4.0 * t + 3.0 * t * t) * cx1 +
+ t * (3.0 * t - 2.0) * cx2) / ds
+
+
+def d2val(ss: float, X: np.ndarray, XS: np.ndarray, S: np.ndarray) -> float:
+ """
+ Estimate the second derivative at ss using central finite differences.
+
+ Parameters
+ ----------
+ ss : float
+ Parameter value.
+ X : ndarray
+ Data array.
+ XS : ndarray
+ Spline derivative array.
+ S : ndarray
+ Parameter array.
+
+ Returns
+ -------
+ float
+ Approximated second derivative.
+ """
+ eps = 1e-6 * (S[-1] - S[0])
+ return (deval(ss + eps, X, XS, S) - deval(ss - eps, X, XS, S)) / (2 * eps)
+
+
+def lefind(X: np.ndarray, XP: np.ndarray, Y: np.ndarray, YP: np.ndarray, S: np.ndarray) -> float:
+ """
+ Locate the leading edge (LE) parameter S_LE. The LE is defined so that the
+ surface tangent is normal to the chord connecting the LE and the trailing edge.
+
+ Parameters
+ ----------
+ X, Y : ndarray
+ Coordinate arrays.
+ XP, YP : ndarray
+ Spline derivative arrays for X and Y.
+ S : ndarray
+ Parameter (arc length) array.
+
+ Returns
+ -------
+ float
+ The parameter value at the leading edge.
+ """
+ N = len(S)
+ dseps = (S[-1] - S[0]) * 1.0e-5
+ x_te = 0.5 * (X[0] + X[-1])
+ y_te = 0.5 * (Y[0] + Y[-1])
+ s_le = None
+ for i in range(2, N - 2):
+ dx_te = X[i] - x_te
+ dy_te = Y[i] - y_te
+ dx = X[i + 1] - X[i]
+ dy = Y[i + 1] - Y[i]
+ dotp = dx_te * dx + dy_te * dy
+ if dotp < 0.0:
+ s_le = S[i]
+ break
+ if s_le is None:
+ s_le = S[N // 2]
+ i_guess = np.searchsorted(S, s_le)
+ if i_guess > 0 and np.isclose(S[i_guess], S[i_guess - 1]):
+ return S[i_guess]
+ for _ in range(50):
+ x_le = seval(s_le, X, XP, S)
+ y_le = seval(s_le, Y, YP, S)
+ dxds = deval(s_le, X, XP, S)
+ dyds = deval(s_le, Y, YP, S)
+ dxdd = d2val(s_le, X, XP, S)
+ dydd = d2val(s_le, Y, YP, S)
+ x_chord = x_le - x_te
+ y_chord = y_le - y_te
+ res = x_chord * dxds + y_chord * dyds
+ ress = dxds**2 + dyds**2 + x_chord * dxdd + y_chord * dydd
+ ds_le = -res / ress
+ ds_le = max(ds_le, -0.02 * abs(x_chord + y_chord))
+ ds_le = min(ds_le, 0.02 * abs(x_chord + y_chord))
+ s_le += ds_le
+ if abs(ds_le) < dseps:
+ return s_le
+ print("LEFIND: LE point not found. Continuing...")
+ return s_le
+
+
+def sopps(
+ si: float,
+ X: np.ndarray,
+ XP: np.ndarray,
+ Y: np.ndarray,
+ YP: np.ndarray,
+ S: np.ndarray,
+ s_le: float
+) -> float:
+ """
+ Find the opposite point parameter for a given point at parameter si.
+ The opposite point is determined by matching the chordwise coordinate.
+
+ Parameters
+ ----------
+ si : float
+ Parameter value of the given point.
+ X, Y : ndarray
+ Coordinate arrays.
+ XP, YP : ndarray
+ Spline derivative arrays.
+ S : ndarray
+ Parameter (arc-length) array.
+ s_le : float
+ Leading edge parameter.
+
+ Returns
+ -------
+ float
+ Parameter of the opposite point.
+ """
+ N = len(S)
+ s_len = S[-1] - S[0]
+ x_le = seval(s_le, X, XP, S)
+ y_le = seval(s_le, Y, YP, S)
+ x_te = 0.5 * (X[0] + X[-1])
+ y_te = 0.5 * (Y[0] + Y[-1])
+ chord = np.sqrt((x_te - x_le) ** 2 + (y_te - y_le) ** 2)
+ dxc = (x_te - x_le) / chord
+ dyc = (y_te - y_le) / chord
+ if si < s_le:
+ i_idx = 0
+ i_opp = N - 1
+ else:
+ i_idx = N - 1
+ i_opp = 0
+ denom = S[i_idx] - s_le
+ s_frac = (si - s_le) / denom if denom != 0 else 0.0
+ s_opp = s_le + s_frac * (S[i_opp] - s_le)
+ if abs(s_frac) <= 1.0e-5:
+ return s_le
+ xi = seval(si, X, XP, S)
+ yi = seval(si, Y, YP, S)
+ x_le = seval(s_le, X, XP, S)
+ y_le = seval(s_le, Y, YP, S)
+ xbar = (xi - x_le) * dxc + (yi - y_le) * dyc
+ for _ in range(12):
+ x_opp = seval(s_opp, X, XP, S)
+ y_opp = seval(s_opp, Y, YP, S)
+ xopp_d = deval(s_opp, X, XP, S)
+ yopp_d = deval(s_opp, Y, YP, S)
+ res = (x_opp - x_le) * dxc + (y_opp - y_le) * dyc - xbar
+ resd = xopp_d * dxc + yopp_d * dyc
+ if abs(res) / s_len < 1.0e-5:
+ break
+ if resd == 0.0:
+ print("SOPPS: Opposite-point location failed. Continuing...")
+ s_opp = s_le + s_frac * (S[i_opp] - s_le)
+ break
+ ds_opp = -res / resd
+ s_opp += ds_opp
+ if abs(ds_opp) / s_len < 1.0e-5:
+ break
+ return s_opp
+
+
+def scalc(X: np.ndarray, Y: np.ndarray) -> np.ndarray:
+ """
+ Compute the cumulative arc-length array for the 2-D curve defined by (X, Y).
+
+ Parameters
+ ----------
+ X, Y : ndarray
+ Coordinate arrays.
+
+ Returns
+ -------
+ ndarray
+ Cumulative arc-length.
+ """
+ N = len(X)
+ S = np.empty(N)
+ S[0] = 0.0
+ for i in range(1, N):
+ S[i] = S[i - 1] + np.hypot(X[i] - X[i - 1], Y[i] - Y[i - 1])
+ return S
+
+
+def curv(ss: float, X: np.ndarray, XS: np.ndarray, Y: np.ndarray, YS: np.ndarray, S: np.ndarray) -> float:
+ """
+ Compute the curvature of a splined 2-D curve at parameter ss.
+
+ Parameters
+ ----------
+ ss : float
+ Parameter value.
+ X, Y : ndarray
+ Coordinate arrays.
+ XS, YS : ndarray
+ Spline derivative arrays.
+ S : ndarray
+ Parameter (arc-length) array.
+
+ Returns
+ -------
+ float
+ Curvature.
+ """
+ N = len(S)
+ i = np.searchsorted(S, ss, side='right')
+ if i == 0:
+ i = 1
+ if i >= N:
+ i = N - 1
+ ds = S[i] - S[i - 1]
+ t = (ss - S[i - 1]) / ds
+ f1 = ds * XS[i - 1]
+ f2 = -ds * (2.0 * XS[i - 1] + XS[i]) + 3.0 * (X[i] - X[i - 1])
+ f3 = ds * (XS[i - 1] + XS[i]) - 2.0 * (X[i] - X[i - 1])
+ xd = f1 + t * (2.0 * f2 + t * 3.0 * f3)
+ xdd = 2.0 * f2 + t * 6.0 * f3
+
+ g1 = ds * YS[i - 1]
+ g2 = -ds * (2.0 * YS[i - 1] + YS[i]) + 3.0 * (Y[i] - Y[i - 1])
+ g3 = ds * (YS[i - 1] + YS[i]) - 2.0 * (Y[i] - Y[i - 1])
+ yd = g1 + t * (2.0 * g2 + t * 3.0 * g3)
+ ydd = 2.0 * g2 + t * 6.0 * g3
+
+ return (xd * ydd - yd * xdd) / np.sqrt((xd**2 + yd**2) ** 3)
+
+
+def ae_calc(N: int, X: np.ndarray, Y: np.ndarray, T: np.ndarray, itype: int):
+ """
+ Calculate geometric properties of the shape defined by X, Y.
+
+ Parameters
+ ----------
+ N : int
+ Number of points.
+ X, Y : ndarray
+ Coordinate arrays.
+ T : ndarray
+ Skin-thickness array (used only if itype == 2).
+ itype : int
+ 1 to integrate over whole area (dx dy); 2 to integrate over skin area (T ds).
+
+ Returns
+ -------
+ tuple
+ (area, xcen, ycen, ei11, ei22, apx1, apx2)
+ """
+ PI = np.pi
+ sint = 0.0
+ aint = 0.0
+ xint = 0.0
+ yint = 0.0
+ xxint = 0.0
+ xyint = 0.0
+ yyint = 0.0
+ for i in range(N):
+ ip = i + 1 if i < N - 1 else 0
+ dx = X[i] - X[ip]
+ dy = Y[i] - Y[ip]
+ xa = 0.5 * (X[i] + X[ip])
+ ya = 0.5 * (Y[i] + Y[ip])
+ ta = 0.5 * (T[i] + T[ip])
+ ds = np.hypot(dx, dy)
+ sint += ds
+ if itype == 1:
+ da = ya * dx
+ aint += da
+ xint += xa * da
+ yint += ya * da / 2.0
+ xxint += xa**2 * da
+ xyint += xa * ya * da / 2.0
+ yyint += ya**2 * da / 3.0
+ else:
+ da = ta * ds
+ aint += da
+ xint += xa * da
+ yint += ya * da
+ xxint += xa**2 * da
+ xyint += xa * ya * da
+ yyint += ya**2 * da
+ if aint == 0.0:
+ return (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
+ xcen = xint / aint
+ ycen = yint / aint
+ eiixx = yyint - ycen**2 * aint
+ eixy = xyint - xcen * ycen * aint
+ eiiyy = xxint - xcen**2 * aint
+ eisq = 0.25 * (eiixx - eiiyy)**2 + eixy**2
+ sgn = 1.0 if (eiiyy - eiixx) >= 0.0 else -1.0
+ ei11 = 0.5 * (eiixx + eiiyy) - sgn * np.sqrt(eisq)
+ ei22 = 0.5 * (eiixx + eiiyy) + sgn * np.sqrt(eisq)
+ if ei11 == 0.0 or ei22 == 0.0:
+ apx1 = 0.0
+ apx2 = np.arctan2(1.0, 0.0)
+ elif eisq / (ei11 * ei22) < (0.001 * sint)**4:
+ apx1 = 0.0
+ apx2 = np.arctan2(1.0, 0.0)
+ else:
+ c1 = eixy
+ s1 = eiixx - ei11
+ c2 = eixy
+ s2 = eiixx - ei22
+ if abs(s1) > abs(s2):
+ apx1 = np.arctan2(s1, c1)
+ apx2 = apx1 + 0.5 * PI
+ else:
+ apx2 = np.arctan2(s2, c2)
+ apx1 = apx2 - 0.5 * PI
+ if apx1 < -0.5 * PI:
+ apx1 += PI
+ if apx1 > 0.5 * PI:
+ apx1 -= PI
+ if apx2 < -0.5 * PI:
+ apx2 += PI
+ if apx2 > 0.5 * PI:
+ apx2 -= PI
+ return aint, xcen, ycen, ei11, ei22, apx1, apx2
+
+
+def tc_calc(X: np.ndarray, XP: np.ndarray, Y: np.ndarray, YP: np.ndarray, S: np.ndarray):
+ """
+ Compute the maximum thickness and camber from discrete airfoil points.
+
+ Parameters
+ ----------
+ X, Y : ndarray
+ Coordinate arrays.
+ XP, YP : ndarray
+ Spline derivative arrays.
+ S : ndarray
+ Parameter (arc-length) array.
+
+ Returns
+ -------
+ tuple
+ (thickness, x_at_thickness, camber, x_at_camber)
+ """
+ s_le = lefind(X, XP, Y, YP, S)
+ x_le = seval(s_le, X, XP, S)
+ y_le = seval(s_le, Y, YP, S)
+ x_te = 0.5 * (X[0] + X[-1])
+ y_te = 0.5 * (Y[0] + Y[-1])
+ chord = np.hypot(x_te - x_le, y_te - y_le)
+ dxc = (x_te - x_le) / chord
+ dyc = (y_te - y_le) / chord
+ thickness = 0.0
+ x_thick = 0.0
+ camber = 0.0
+ x_camber = 0.0
+ N = len(S)
+ for i in range(N):
+ x_bar = (X[i] - x_le) * dxc + (Y[i] - y_le) * dyc
+ y_bar = (Y[i] - y_le) * dxc - (X[i] - x_le) * dyc
+ s_opp = sopps(S[i], X, XP, Y, YP, S, s_le)
+ x_opp = seval(s_opp, X, XP, S)
+ y_opp = seval(s_opp, Y, YP, S)
+ y_bar_opp = (y_opp - y_le) * dxc - (x_opp - x_le) * dyc
+ yc = 0.5 * (y_bar + y_bar_opp)
+ yt = abs(y_bar - y_bar_opp)
+ if abs(yc) > abs(camber):
+ camber = yc
+ x_camber = x_opp
+ if abs(yt) > abs(thickness):
+ thickness = yt
+ x_thick = x_opp
+ return thickness, x_thick, camber, x_camber
+
+
+def geopar(
+ X: np.ndarray, XP: np.ndarray, Y: np.ndarray, YP: np.ndarray,
+ S: np.ndarray, T: np.ndarray, verbose: bool=False
+) -> dict:
+ """
+ Compute overall geometric parameters of the airfoil shape.
+
+ Parameters
+ ----------
+ X, Y : ndarray
+ Coordinate arrays.
+ XP, YP : ndarray
+ Spline derivative arrays.
+ S : ndarray
+ Parameter (arc-length) array.
+ T : ndarray
+ Skin-thickness array (usually ones).
+
+ Returns
+ -------
+ dict
+ Dictionary of geometric parameters.
+ """
+ N = len(S)
+ s_le = lefind(X, XP, Y, YP, S)
+ x_le = seval(s_le, X, XP, S)
+ y_le = seval(s_le, Y, YP, S)
+ x_te = 0.5 * (X[0] + X[-1])
+ y_te = 0.5 * (Y[0] + Y[-1])
+ chord = np.hypot(x_te - x_le, y_te - y_le)
+ curv_le = curv(s_le, X, XP, Y, YP, S)
+ rad_le = 1.0 / curv_le if abs(curv_le) > 0.001 * (S[-1] - S[0]) else 0.0
+ ang1 = np.arctan2(-YP[0], -XP[0])
+ ang2 = np.arctan2(YP[-1], XP[-1])
+ ang_te = ang2 - ang1
+ area, xcena, ycena, ei11a, ei22a, apx1a, apx2a = ae_calc(N, X, Y, T, itype=1)
+ slen, xcent, ycent, ei11t, ei22t, apx1t, apx2t = ae_calc(N, X, Y, T, itype=2)
+ thick, x_thick, cambr, x_camber = tc_calc(X, XP, Y, YP, S)
+ if verbose:
+ print(f"Max thickness = {thick: .6f} at x = {x_thick: .3f}")
+ print(f"Max camber = {cambr: .6f} at x = {x_camber: .3f}")
+ return {
+ "SLE": s_le, "CHORD": chord, "AREA": area, "RADLE": rad_le, "ANGTE": ang_te,
+ "EI11A": ei11a, "EI22A": ei22a, "APX1A": apx1a, "APX2A": apx2a,
+ "EI11T": ei11t, "EI22T": ei22t, "APX1T": apx1t, "APX2T": apx2t,
+ "THICK": thick, "CAMBR": cambr
+ }
+
+
+def sortol(x_arr: np.ndarray, y_arr: np.ndarray, tol: float) -> (np.ndarray, np.ndarray):
+ """
+ Sort (x, y) pairs by x.
+
+ Parameters
+ ----------
+ x_arr : ndarray
+ x-coordinate array.
+ y_arr : ndarray
+ y-coordinate array.
+ tol : float
+ Tolerance (unused here, kept for compatibility).
+
+ Returns
+ -------
+ tuple
+ Sorted (x_arr, y_arr) arrays.
+ """
+ order = np.argsort(x_arr)
+ return x_arr[order], y_arr[order]
+
+
+def getcam(
+ X: np.ndarray, XP: np.ndarray, Y: np.ndarray, YP: np.ndarray, S: np.ndarray
+) -> (np.ndarray, np.ndarray, int, np.ndarray, np.ndarray, int):
+ """
+ Compute the camber line and thickness distribution for an airfoil.
+
+ Parameters
+ ----------
+ X, Y : ndarray
+ Coordinate arrays.
+ XP, YP : ndarray
+ Spline derivative arrays.
+ S : ndarray
+ Parameter (arc-length) array.
+
+ Returns
+ -------
+ tuple
+ (x_camber, y_camber, n_camber, x_thick, y_thick, n_thick)
+ """
+ N = len(S)
+ s_le = lefind(X, XP, Y, YP, S)
+ x_le = seval(s_le, X, XP, S)
+ y_le = seval(s_le, Y, YP, S)
+ x_camber = np.empty(N)
+ y_camber = np.empty(N)
+ x_thick = np.empty(N)
+ y_thick = np.empty(N)
+ for i in range(N):
+ if i == 0:
+ x_opp = X[-1]
+ y_opp = Y[-1]
+ elif i == N - 1:
+ x_opp = X[0]
+ y_opp = Y[0]
+ else:
+ s_opp = sopps(S[i], X, XP, Y, YP, S, s_le)
+ x_opp = seval(s_opp, X, XP, S)
+ y_opp = seval(s_opp, Y, YP, S)
+ x_camber[i] = 0.5 * (X[i] + x_opp)
+ y_camber[i] = 0.5 * (Y[i] + y_opp)
+ x_thick[i] = 0.5 * (X[i] + x_opp)
+ y_thick[i] = abs(0.5 * (Y[i] - y_opp))
+ # Append the leading-edge point
+ x_camber = np.append(x_camber, x_le)
+ y_camber = np.append(y_camber, y_le)
+ n_camber = len(x_camber)
+ x_thick = np.append(x_thick, x_le)
+ y_thick = np.append(y_thick, 0.0)
+ n_thick = len(x_thick)
+ tol_val = 1.0e-5 * (S[-1] - S[0])
+ x_camber, y_camber = sortol(x_camber, y_camber, tol_val)
+ y_camber -= y_camber[0]
+ x_thick, y_thick = sortol(x_thick, y_thick, tol_val)
+ return x_camber, y_camber, n_camber, x_thick, y_thick, n_thick
+
+
+def thkcam(
+ tfac: float,
+ cfac: float,
+ XB: np.ndarray,
+ XBP: np.ndarray,
+ YB: np.ndarray,
+ YBP: np.ndarray,
+ SB: np.ndarray
+) -> (np.ndarray, np.ndarray, dict):
+ """
+ Adjust the airfoil (buffer) thickness and camber.
+
+ Parameters
+ ----------
+ tfac : float
+ Thickness scaling factor.
+ cfac : float
+ Camber scaling factor.
+ XB, YB : ndarray
+ Buffer airfoil coordinate arrays.
+ XBP, YBP : ndarray
+ Buffer airfoil spline derivative arrays.
+ SB : ndarray
+ Parameter (arc-length) array for the buffer.
+
+ Returns
+ -------
+ tuple
+ Updated XB, YB arrays and a dictionary of geometric parameters.
+ """
+ NB = len(SB)
+ s_le = lefind(XB, XBP, YB, YBP, SB)
+ x_le = seval(s_le, XB, XBP, SB)
+ y_le = seval(s_le, YB, YBP, SB)
+ x_te = 0.5 * (XB[0] + XB[-1])
+ y_te = 0.5 * (YB[0] + YB[-1])
+ chord = np.hypot(x_te - x_le, y_te - y_le)
+ dxc = (x_te - x_le) / chord
+ dyc = (y_te - y_le) / chord
+ W1 = np.empty(NB)
+ W2 = np.empty(NB)
+ for i in range(NB):
+ if i == 0:
+ xb_opp = XB[-1]
+ yb_opp = YB[-1]
+ elif i == NB - 1:
+ xb_opp = XB[0]
+ yb_opp = YB[0]
+ else:
+ s_opp = sopps(SB[i], XB, XBP, YB, YBP, SB, s_le)
+ xb_opp = seval(s_opp, XB, XBP, SB)
+ yb_opp = seval(s_opp, YB, YBP, SB)
+ xc_avg = 0.5 * (XB[i] + xb_opp) * dxc + 0.5 * (YB[i] + yb_opp) * dyc
+ yc_avg = cfac * (0.5 * (YB[i] + yb_opp) * dxc - 0.5 * (XB[i] + xb_opp) * dyc)
+ xc_del = 0.5 * (XB[i] - xb_opp) * dxc + 0.5 * (YB[i] - yb_opp) * dyc
+ yc_del = tfac * (0.5 * (YB[i] - yb_opp) * dxc - 0.5 * (XB[i] - xb_opp) * dyc)
+ W1[i] = (xc_avg + xc_del) * dxc - (yc_avg + yc_del) * dyc
+ W2[i] = (yc_avg + yc_del) * dxc + (xc_avg + xc_del) * dyc
+
+ # Update buffer coordinates with new values
+ XB[:] = W1
+ YB[:] = W2
+
+ # Recalculate arc-length and update spline derivatives
+ SB_new = scalc(XB, YB)
+ XBP = segspl(XB, SB_new)
+ YBP = segspl(YB, SB_new)
+ # For GEOPAR, define a thickness array (using ones)
+ T = np.ones_like(XB)
+ geo_params = geopar(XB, XBP, YB, YBP, SB_new, T)
+ return XB, YB, geo_params
+
+def load_profile(file_path: str=""):
+ """
+ Reads an airfoil profile file and returns the x and y coordinates.
+ Assumes the file contains two whitespace-separated columns.
+ """
+ try:
+ data = np.loadtxt(file_path)
+ except Exception as e:
+ print(f"Error reading {file_path}: {e}")
+ return None, None
+ if data.ndim == 1 or data.shape[1] < 2:
+ print(f"File {file_path} does not have two columns.")
+ return None, None
+
+ return data # x, y coordinates
+
+
+def rescale_profile(coordinates: np.ndarray=None, scale_t_to_c: float=1.0, scale_camber: float=1.0) -> np.ndarray:
+ """ Rescale profile according to thickness to chord or to camber
+
+ Keyword Arguments:
+ coordinates {np.ndarray} -- coordinates x, y (from TE,UP -> LE -> TE,DN) (default: {None})
+ scale_t_to_c {float} -- scaling factor to scale current thickness to chord (default: {1.0})
+ scale_camber {float} -- scaling factor to scale current camber (default: {1.0})
+
+ Returns:
+ np.ndarray -- scaled coordinates
+ """
+ x = coordinates[:, 0]
+ y = coordinates[:, 1]
+
+ # Spline points
+ s = np.linspace(0, 1, coordinates.shape[0])
+ xp = segspl(x, s)
+ yp = segspl(y, s)
+
+ xcpy = x.copy()
+ ycpy = y.copy()
+ scpy = s.copy()
+ xpcpy = xp.copy()
+ ypcpy = yp.copy()
+ x_new, y_new, _ = thkcam(scale_t_to_c, scale_camber, xcpy, xpcpy, ycpy, ypcpy, scpy)
+ return np.column_stack((x_new, y_new))
+
+def save_profile(coordinates: np.ndarray=None, file: str="./scaled_profile.dat") -> None:
+ """ Save profile coordinates in file
+
+ Keyword Arguments:
+ coordinates {np.ndarray} -- profile coordinates (default: {None})
+ file {str} -- file + (default: {"./scaled_profile.dat"})
+ """
+ np.savetxt(file, coordinates, fmt="%.8f", delimiter=" ")
+ return
\ No newline at end of file
diff --git a/pyavlunicadopackage/pyproject.toml b/pyavlunicadopackage/pyproject.toml
new file mode 100644
index 0000000000000000000000000000000000000000..ab380c63ed8f549170078113797d7be030c1bbf9
--- /dev/null
+++ b/pyavlunicadopackage/pyproject.toml
@@ -0,0 +1,27 @@
+[build-system]
+requires = ["setuptools>=42", "wheel"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "pyavlunicado"
+version = "1.0.0"
+description = "Unicado avl interface to build generate aerodynamic derivatives."
+authors = [{ name = "Christopher Ruwisch", email = "christopher.ruwisch@tu-berlin.de" }]
+readme = "README.md"
+license = { file = "LICENSE" }
+keywords = ["avl"]
+classifiers = [
+ "Programming Language :: Python :: 3",
+ "License :: OSI Approved :: GNU General Public License v3 or later (GPLv3+)",
+ "Operating System :: OS Independent",
+]
+dependencies = [
+ "numpy>=1.26.4",
+]
+
+[tool.setuptools]
+packages = ["pyavlunicado"]
+
+[project.urls]
+homepage = "https://unicado.pages.rwth-aachen.de/unicado.gitlab.io/"
+repository = "https://git.rwth-aachen.de/unicado/libraries"
\ No newline at end of file