Added script for VSPAERO .history file -> JSBSim FDM table conversion, added sinusoidal method to utils.interpolator

aircraft/vsphist2jsbtable.py

@@ -0,0 +1,123 @@
+#!/usr/bin/env python
+#-*- coding:utf-8 -*-
+
+import argparse
+import os, sys
+from fgtools.utils.interpolator import Interpolator
+
+class Case:
+	def __init__(self, textcase):
+		textcase.pop(0) # get rid of the Name Value Unit legend line
+		self.Sref = float(textcase.pop(0).split()[1])
+		self.Bref = float(textcase.pop(0).split()[1])
+		self.Cref = float(textcase.pop(0).split()[1])
+		self.Xcg = float(textcase.pop(0).split()[1])
+		self.Ycg = float(textcase.pop(0).split()[1])
+		self.Zcg = float(textcase.pop(0).split()[1])
+		self.Mach = float(textcase.pop(0).split()[1])
+		self.AoA = float(textcase.pop(0).split()[1])
+		self.Beta = float(textcase.pop(0).split()[1])
+		self.Rho = float(textcase.pop(0).split()[1])
+		self.Vinf = float(textcase.pop(0).split()[1])
+		self.RollRate = float(textcase.pop(0).split()[1])
+		self.PitchRate = float(textcase.pop(0).split()[1])
+		self.YawRate = float(textcase.pop(0).split()[1])
+		textcase.pop(0) # get rid of the Solver case line
+		textcase.pop(0) # get rid of another legend line
+		lastiterindex = textcase.index("Skin Friction Drag Break Out:") - 1
+		self.CL, self.CDo, self.CDi, self.CDtot, self.CS, self.LD, self.E, self.CFx, self.CFy, self.CFz, self.CMx, self.CMy, self.CMz, self.CDtrefftz, self.TQS = map(float, textcase.pop(lastiterindex).split()[4:]) # we don't need Mach, AoA, Beta again
+	
+	def __str__(self):
+		return "Case(" + ", ".join(str(k) + " = " + str(v) for k, v in vars(self).items()) + ")"
+	
+	def __repr__(self):
+		return self.__str__()
+
+def get_cases(path):
+	path = os.path.abspath(path)
+	
+	if not path.endswith(".history"):
+		print("The specified file", path, "is not a VSPAERO .history file - exiting")
+		sys.exit(1)
+	
+	with open(path, "r") as histf:
+		textcases = []
+		lines = histf.readlines()[1:]
+		textcase = []
+		cases = {}
+		
+		for line in lines:
+			if line.startswith("****"):
+				textcases.append(list(filter(None, textcase)))
+				textcase = []
+				continue
+			textcase.append(line.strip())
+		textcases.append(list(filter(None, textcase)))
+		
+		for textcase in textcases:
+			case = Case(textcase)
+			if not case.AoA in cases:
+				cases[case.AoA] = {}
+			cases[case.AoA][case.Beta] = case
+	
+	return cases
+
+def print_table(cases, coeff, indent, precision):
+	coeffs = {}
+	for AoA in cases:
+		for Beta in cases[AoA]:
+			if not hasattr(cases[AoA][Beta], coeff):
+				print("The coefficient", coeff, "does not exist in the specified .history file - exiting")
+				sys.exit(1)
+			
+			if not AoA in coeffs:
+				coeffs[AoA] = {}
+			coeffs[AoA][Beta] = getattr(cases[AoA][Beta], coeff)
+	
+	print("<table>")
+	print(indent + '<independentVar lookup="row">aero/alpha-deg</independentVar>')
+	print(indent + '<independentVar lookup="column">aero/beta-deg</independentVar>')
+	#print(indent + '<independentVar lookup="table">velocities/mach</independentVar>')
+	print(indent + "<tableData>")
+	print(indent + indent + indent + indent.join(map(str, coeffs[list(coeffs.keys())[0]].keys())))
+	for AoA in coeffs:
+		print(indent + indent + str(AoA), end="")
+		for Beta in coeffs[AoA]:
+			print(indent + ("%." + str(precision) + "f") % coeffs[AoA][Beta], end="")
+		print()
+	print(indent + "</tableData>")
+	print("</table>")
+
+
+if __name__ == "__main__":
+	argp = argparse.ArgumentParser(description="vsphist2jsbtable.py - Takes a VSPAERO .history file and a coefficient as input and outputs a JSBSim interpolation table for that coefficient from the VSPAERO cases")
+	
+	argp.add_argument(
+		"-c", "--coeff",
+		help="Name of the coefficient to produce a table for",
+		required=True,
+	)
+	
+	argp.add_argument(
+		"--indentation",
+		help="The argument of this option will be used as one level of indentation, defaults to a tab",
+		default="\t"
+	)
+	
+	argp.add_argument(
+		"-p", "--precision",
+		help="How many decimal places the numbers in the table should have, defaults to 6",
+		type=int,
+		default=6
+	)
+	
+	argp.add_argument(
+		"input_file",
+		help="VSPAERO .history file",
+	)
+	
+	args = argp.parse_args()
+	
+	cases = get_cases(args.input_file)
+	print_table(cases, args.coeff, args.indentation, args.precision)
+

utils/interpolator.py

@@ -1,14 +1,14 @@
 #!/usr/bin/env python
 #-*- coding:utf-8 -*-
 
+import math
+
 class Interpolator:
 	def __init__(self):
 		self._indexes = []
 		self._values = []
 		self._sorted = False
 	
-		self.methods = {"linear": self._interpolate_linear}
-	
 	def add_value(self, index, value):
 		if type(index) not in (int, float) or type(value) not in (int, float):
 			try:
@@ -26,7 +26,7 @@ class Interpolator:
 			self.add_value(i, v)
 	
 	def interpolate(self, index, extrapolate=True, method="linear", sort=True):
-		if not method in self.methods:
+		if not hasattr(self, "_interpolate_" + method):
 			raise NotImplementedError(f"Interpolator.interpolate: interpolation method '{method}' not yet supported")
 		
 		if len(self._indexes) < 2:
@@ -38,7 +38,22 @@ class Interpolator:
 			self._values.sort()
 			self._sorted = True
 		
-		return self.methods[method](index, extrapolate)
+		if index in self._indexes:
+			return self._values[self._indexes.index(index)]
+		
+		if self._indexes[0] < index < self._indexes[-1]:
+			return getattr(self, "_interpolate_" + method)(index, extrapolate)
+		else:
+			if not extrapolate:
+				if index < self._indexes[0]:
+					return self._values[0]
+				else:
+					return self._values[-1]
+			else:
+				if index < self._indexes[0]:
+					return self._values[1] + (index - self._indexes[1]) / (self._indexes[0] - self._indexes[1]) * (self._values[0] - self._values[1])
+				else:
+					return self._values[-2] + (index - self._indexes[-2]) / (self._indexes[-1] - self._indexes[-2]) * (self._values[-1] - self._values[-2])
 	
 	def _find_neighbours(self, index):
 		lower = upper = 0
@@ -53,30 +68,24 @@ class Interpolator:
 		return lower, upper
 	
 	def _interpolate_linear(self, index, extrapolate=True):
-		if index in self._indexes:
+		lower, upper = self._find_neighbours(index)
-			return self._values[self._indexes.index(index)]
+		return self._values[lower] + (self._values[upper] - self._values[lower]) * (index - self._indexes[lower]) / (self._indexes[upper] - self._indexes[lower])
 	
-		if self._indexes[0] < index < self._indexes[-1]:
+	def _interpolate_sinusoidal(self, index, extrapolate=True):
-			lower, upper = self._find_neighbours(index)
+		lower, upper = self._find_neighbours(index)
-			return self._values[lower] + (self._values[upper] - self._values[lower]) * (index - self._indexes[lower]) / (self._indexes[upper] - self._indexes[lower]) 
+		return self._values[lower] + (self._values[upper] - self._values[lower]) * math.sin(math.radians((index - self._indexes[lower]) / (self._indexes[upper] - self._indexes[lower]) * 90))
-		else:
-			if not extrapolate:
-				if index < self._indexes[0]:
-					return self._values[0]
-				else:
-					return self._values[-1]
-			else:
-				if index < self._indexes[0]:
-					return self._values[1] + (index - self._indexes[1]) / (self._indexes[0] - self._indexes[1]) * (self._values[0] - self._values[1])
-				else:
-					return self._values[-2] + (index - self._indexes[-2]) / (self._indexes[-1] - self._indexes[-2]) * (self._values[-1] - self._values[-2])
 
 # run test if run directly
 if __name__ == "__main__":
-	print("Test results")
+	print("Test results:")
+	print()
 	i = Interpolator()
 	i.add_values((0, 10, 20), (0, 20, 30))
-	for test_val in (-5, 0, 1, 2, 3.5, 5.55555, 9, 10, 15, 100):
+	test_vals = (-5, 0, 1, 2, 3.5, 5.55555, 10, 15, 35, 100)
-		print(test_val, i.interpolate(test_val))
 	
+	for method in ("linear", "sinusoidal"):
+		print(method.capitalize() + ":")
+		for test_val in test_vals:
+			print(test_val, "=>", i.interpolate(test_val, method=method))
+		print()