Skip to content

Instantly share code, notes, and snippets.

@fredvol

fredvol/share_stu.py

Created December 7, 2023 13:30
Show Gist options
  • Save fredvol/3a6267854bce52acd9ea9de98f05d297 to your computer and use it in GitHub Desktop.
Save fredvol/3a6267854bce52acd9ea9de98f05d297 to your computer and use it in GitHub Desktop.
Download ZIP
model for aerosandbox
Raw
share_stu.py
# MODEL FILE FOR QUICK 3D POLAR CALCULATION
# 17/07/2023 - frd
# %%
import json
import os
import re
from typing import List, Union
import aerosandbox as asb
import aerosandbox.numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from scipy.linalg import expm, norm
from copy import deepcopy
from skspatial.objects import Line
import plotly.graph_objects as go
import plotly.express as px
# %%
glider_model_path = os.path.join(
"ref_files", "Enzo3_M_production_weightCost_state.json"
)
# %% FUNCTIONS
def distance_between_2pts(point1, point2):
# Convert input points to NumPy arrays for vector operations
point1 = np.array(point1)
point2 = np.array(point2)
# Calculate the difference between the two points
difference = point2 - point1
# Calculate the Euclidean distance
distance = np.linalg.norm(difference)
return distance
def convert_float_list(mstr):
"""
Convert a comma-separated string of numbers to a list of floats.
Used to convert a list of coord of 3D point to a python list.
Parameters:
mstr (str): A string containing comma-separated numerical values.
Returns:
list of float: A list containing the numerical values as floats.
Example:
>>> convert_float_list("1.5, 2.3, 3.7")
[1.5, 2.3, 3.7]
Note:
- If the input `mstr` is not a string, the function will raise a TypeError.
- The function does not perform extensive error-checking on the input string.
- The input string must contain numerical values separated by commas, with no spaces.
- Any invalid or improperly formatted values may result in unexpected behavior.
"""
if type(mstr) is str:
return [float(i) for i in mstr.split(",")]
def unit_vector(vector):
"""Returns the unit vector of the vector."""
return vector / np.linalg.norm(vector)
def vector_between2pts(point1, point2, unitize=True):
# Convert the input points to NumPy arrays
point1 = np.array(point1)
point2 = np.array(point2)
# Calculate the vector between the two points
vector = point2 - point1
if unitize:
vector = unit_vector(vector)
return vector
def translate_point(point, vector):
# Convert the input point and vector to NumPy arrays
point = np.array(point)
vector = np.array(vector)
# Translate the point by adding the vector
translated_point = point + vector
return translated_point
def convert_mm_to_m(my_list):
return [i / 1000 for i in my_list]
def find_point0_TE(original_point, angle_degrees, distance):
"""0nly valid for root chord, TODO use the referenciel of each Xsecs
Args:
original_point (List): _description_
angle_degrees (float): _description_
distance (float): _description_
Returns:
List of float : coorddinate
"""
# Convert input point to a NumPy array
original_point = np.array(original_point)
# Convert angle from degrees to radians
angle_radians = np.radians(angle_degrees)
# Calculate the new point's coordinates
new_x = original_point[0] + distance * np.cos(angle_radians)
new_z = original_point[1] + distance * np.sin(angle_radians)
# Keep the original Z coordinate
new_y = original_point[2]
return [new_x, new_y, new_z]
def angle_between(v1, v2=(0, 1, 0)):
"""Returns the angle in radians between vectors 'v1' and 'v2'::
v2, default is Y axis
>>> angle_between((1, 0, 0), (0, 1, 0))
1.5707963267948966
>>> angle_between((1, 0, 0), (1, 0, 0))
0.0
>>> angle_between((1, 0, 0), (-1, 0, 0))
3.141592653589793
"""
v1_u = unit_vector(v1)
v2_u = unit_vector(v2)
return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))
def rotate_pt(orign_pt, theta_deg, axis=[0, -1, 0]):
"""
Rotate a 3D point around an arbitrary axis.
The function rotates a 3D point `origin_pt` by an angle `theta_deg` (in degrees)
around an arbitrary axis defined by `axis`.
Parameters:
origin_pt (list or numpy.array): A list or numpy array representing the 3D point to rotate
in the format [x, y, z].
theta_deg (float): The angle of rotation in degrees.
axis (list or numpy.array, optional): A list or numpy array representing the axis of rotation.
The default axis is Y , [0, -1, 0].
Returns:
numpy.array: A numpy array representing the rotated 3D point.
Example:
>>> rotate_pt(point=[1.0, 0.0, 0.0] , angle= 90,axis=[0, 0, 1] )
[0 , 1 , 0]
Note:
- The `axis` should be a unit vector; otherwise, the function will normalize it.
- The function uses the Rodrigues' rotation formula to perform the rotation.
"""
theta = np.radians(theta_deg)
M0 = expm(np.cross(np.eye(3), axis / norm(axis) * theta))
return np.dot(M0, orign_pt)
def convert_sentence_to_dict(sentence, multiplicator_factor=1):
"""
Convert a sentence with key-value pairs into a Python dictionary ( generated chatGPT).
The function takes a sentence containing key-value pairs and converts it into a Python
dictionary where keys are strings and values are lists of three floating-point numbers.
The sentence must have the following format for each key-value pair:
'key:[x, y, z]', where 'key' is a string, and 'x', 'y', 'z' are floating-point numbers.
Parameters:
sentence (str): A sentence containing key-value pairs in the specified format.
multiplicator_factor (float) : value to scale de data from mm to m.
Returns:
dict: A Python dictionary where keys are strings and values are lists of three floating-point numbers.
Note:
- The function assumes the provided sentence follows the specified format for each key-value pair.
- The function does not perform extensive error-checking on the input sentence.
- Any invalid or improperly formatted key-value pairs may result in unexpected behavior.
"""
# Regular expression pattern to match key-value pairs
pattern = r"(\w+):\[(\-*\d+\.\d+),(\-*\d+\.\d+),(\-*\d+\.\d+)\]"
# Find all key-value pairs in the sentence
matches = re.findall(pattern, sentence)
# Create the dictionary
result_dict = {
key: [
float(x) * multiplicator_factor,
float(y) * multiplicator_factor,
float(z) * multiplicator_factor,
]
for key, x, y, z in matches
}
return result_dict
def extract_profile(s):
"""convert profile coord as list space seaprated to numpy N2x array.
the profile is scaled to 0,1 coord.
Args:
s (str): string input
Returns:
np.ndarray: Scaled profile coords
Example:
>>> extract_profile("2068.0,0.0 2034.9,5.7 1971.6,16.7 1895.3,29.9 1813.0,44.1")
array([[1. , 0. ],
[0.9839942 , 0.00275629],
[0.95338491, 0.00807544],
[0.91648936, 0.01445841],
[0.87669246, 0.02132495]])
"""
x = []
y = []
for ligne in s.split(" "):
xi, yi = ligne.split(",")
x.append(float(xi))
y.append(float(yi))
ratio_scale = 1 / max(x)
xs = np.array(x) * ratio_scale
ys = np.array(y) * ratio_scale
return np.array(list(zip(xs, ys)))
def weighted_barycentric_average(points, weights):
"""
Calculate the weighted barycentric average of a list of 3D points. ( generated chatGPT)
The function computes the weighted barycentric average of the given 3D points
using the provided weights.
Parameters:
points (list of lists): A list containing 3D points represented as lists of three
floating-point numbers [x, y, z].
weights (list of float): A list of weights corresponding to each point in the `points` list.
Returns:
list: A list containing the x, y, and z coordinates of the weighted barycentric average.
Raises:
ValueError: If the number of points and weights are not the same
Example:
>>> points = [
... [83.4, 320.1, -7871.3],
... [83.4, 320.1, -7871.3],
... [83.4, 320.1, -7871.3],
... [627.8, 5017.5, -3097.9],
... [627.8, 5017.5, -3097.9],
... ]
>>> weights = [0.009547, 0.006893, 0.009810, 0.000106, 0.000109]
>>> result = weighted_barycentric_average(points, weights)
>>> print(result)
[83.40001397960151, 320.1000336941283, -7871.299988324128]
Note:
The function assumes that the number of points and weights is the same,
"""
if len(points) != len(weights):
raise ValueError("Number of points and weights must be the same.")
weighted_sum = [0.0, 0.0, 0.0]
total_weight = 0.0
for i, point in enumerate(points):
weight = weights[i]
total_weight += weight
for j in range(3):
weighted_sum[j] += point[j] * weight
barycentric_average = [coord_sum / total_weight for coord_sum in weighted_sum]
return barycentric_average
class ParagliderData:
def __init__(self, glider_model_path):
self.source_path = glider_model_path
self.glider_mass = 5 # kg
with open(self.source_path) as f:
# data = json.load(f) # ,object_hook=Generic.from_dict)
# Parse JSON into an object with attributes corresponding to dict keys.
self.original_data = json.loads(
f.read().replace(
"/*Already listed*/", ""
) # , object_hook=lambda d: SimpleNamespace(**d)
)
self.compute_general()
self.compute_ribs()
self.compute_lines()
print(" end importing data from Json file")
def compute_general(self):
# % Import general
interesting_fields = [
"Cells",
"FlatSpan",
"HarnessDragArea",
"Name",
"RiserLength",
"RiserWebs",
"RootChord",
"RootThickness",
"ShortName",
"Size",
"Span",
"VertDist",
"Area",
"HasCentreCell",
"HasCentreRib",
"AspectRatio",
"PilotLoad",
"PilotVerticalCentreOfGravity",
"HarnessDragArea",
]
general_data = dict(
(k, self.original_data[k])
for k in interesting_fields
if k in self.original_data
)
general_data["RisersDict"] = convert_sentence_to_dict(
general_data["RiserWebs"], multiplicator_factor=0.001
)
# convert to numeric
for c in ["PilotLoad", "HarnessDragArea"]:
general_data[c] = float(general_data[c])
# convert mm to m
for c in [
"FlatSpan",
"RiserLength",
"RootChord",
"Span",
"VertDist",
"PilotVerticalCentreOfGravity",
]:
general_data[c] = float(general_data[c]) / 1000
self.general = general_data
print(" General data computed")
def compute_ribs(self):
if hasattr(self, "general"):
# % WORK ON RIBS
ribs = [
item["Children"]
for item in self.original_data["Children"]
if "FoilRibs" in item["$Type"]
][0]
df_ribs_full = pd.json_normalize(ribs)
df_ribs = df_ribs_full[
[
"Name",
"ChordAngle",
"Chord",
"LEPoint",
"TEPoint",
"Profile2D",
"Tabs",
]
]
df_ribs.loc[:, "ChordAngle"] = pd.to_numeric(df_ribs["ChordAngle"])
df_ribs.loc[:, "Chord"] = pd.to_numeric(df_ribs["Chord"])
df_ribs.loc[:, "Chord"] /= 1000 # mm to m conversion
df_ribs.loc[:, "LEPoint"] = df_ribs["LEPoint"].apply(convert_float_list)
df_ribs.loc[:, "LEPoint"] = df_ribs["LEPoint"].apply(convert_mm_to_m)
df_ribs.loc[:, "TEPoint"] = df_ribs["TEPoint"].apply(convert_float_list)
df_ribs.loc[:, "TEPoint"] = df_ribs["TEPoint"].apply(convert_mm_to_m)
if self.general["HasCentreCell"]:
# add a virtual center rib
Rib0 = deepcopy(df_ribs.loc[0])
le_point0 = list(Rib0["LEPoint"])
le_point0[0] = 0
le_point0[1] = 0
te_point0 = list(Rib0["TEPoint"])
te_point0[0] = Rib0["TEPoint"][0]
te_point0[1] = 0
Rib0["TEPoint"] = te_point0
df_ribs.loc[-1] = Rib0
df_ribs.index = df_ribs.index + 1
df_ribs = df_ribs.sort_index()
# create local vector
# TODO comapre to Grashopper
# % simplify tabs
def extract_tabs(tabs):
tabs_return = []
if len(tabs) > 0:
for tab_dict in tabs:
if tab_dict["LoadBearing"]:
tabs_return.append(
dict(
(k, tab_dict[k])
for k in ["Type", "PcntChord", "PositionBottom3D"]
if k in tab_dict
)
)
return tabs_return
df_ribs.loc[:, "Tabs"] = df_ribs["Tabs"].apply(extract_tabs)
self.ribs_df = df_ribs
print(" Lines data computed")
else:
print(" impossible general data not computed")
def compute_lines(self):
# % WORK ON Lines
self.line_splice_ratio = 0.5
lines = [
item["Children"]
for item in self.original_data["Children"]
if "SuspensionLines" in item["$Type"]
][0]
df_lines_full = pd.json_normalize(lines)
df_lines_na_include = df_lines_full[
[
"Name",
"Diameter",
"AdjustedLength",
"Material",
"Riser",
"@Point0",
"@Point1",
]
].copy() # here the SettingWithCopyWarning could be triger
df_lines = df_lines_na_include.dropna()
df_lines.loc[:, "Diameter"] = pd.to_numeric(df_lines["Diameter"])
df_lines.loc[:, "AdjustedLength"] = pd.to_numeric(df_lines["AdjustedLength"])
df_lines.loc[:, "AdjustedLength"] /= 1000
df_lines.loc[:, "SpliceLength"] = 0.100 * 2
df_lines.loc[:, "cx"] = 1
df_lines.loc[:, "Scx"] = (
df_lines["Diameter"] / 1000 * df_lines["cx"] * df_lines["AdjustedLength"]
)
df_lines.loc[:, "Scx_splice"] = (
df_lines["Diameter"]
/ 1000
* df_lines["cx"]
* self.line_splice_ratio
* df_lines["SpliceLength"]
)
df_lines.loc[:, "Scx_tot"] = df_lines["Scx_splice"] + df_lines["Scx"]
total_line_drag = df_lines["Scx_tot"].sum() * 2
# compute barycenter
df_lines.loc[:, "@Point0"] = df_lines["@Point0"].apply(convert_float_list)
df_lines.loc[:, "@Point1"] = df_lines["@Point1"].apply(convert_float_list)
df_lines.loc[:, "@Point0"] = df_lines["@Point0"].apply(convert_mm_to_m)
df_lines.loc[:, "@Point1"] = df_lines["@Point1"].apply(convert_mm_to_m)
line_application_point_one_side = weighted_barycentric_average(
list(df_lines["@Point0"]), list(df_lines["Scx_tot"])
)
line_application_point = line_application_point_one_side
line_application_point[1] = 0 # remove y component
self.lines_scx = total_line_drag
self.lines_line_application_point = line_application_point
self.lines_df = df_lines
def generate_wing(self):
# % Create Wing
mainwing_xsecs = []
for i, row in self.ribs_df.iterrows():
xsec = asb.WingXSec(
xyz_le=row["LEPoint"],
chord=row["Chord"],
twist=row["ChordAngle"],
airfoil=asb.Airfoil(
row["Name"], coordinates=extract_profile(row["Profile2D"])
),
)
mainwing_xsecs.append(xsec)
mainwing = asb.Wing(
name="Main",
xsecs=mainwing_xsecs,
symmetric=True,
is_main_wing=True,
is_horizontal_stabilizer=False,
is_vertical_stabilizer=False,
is_fuselage=False,
)
return mainwing
# %%
t = ParagliderData(glider_model_path)
# %%
class Paraglider(asb.Airplane):
def __init__(self, glider_model_path, derotate=True):
# PArse json file
self.data = ParagliderData(glider_model_path)
# Call the constructor of the base class 'aerosandbox.Airplane'
super().__init__(
name=self.data.general["ShortName"],
xyz_ref=[0, 0, 0],
wings=[self.data.generate_wing()],
)
self.data.general["original_angle"] = self.main_ChordAngle()
self.data.general["hasbeenDerotated"] = False
if derotate and self.main_ChordAngle() != 0:
self.derotate_all()
self._compute_main_pts()
# self.data.[elt.twist for elt in glider.wings[0].xsecs]
self._store_twist_origin()
def draw_skeleton(self):
df = self.data.ribs_df
df[["X_LE", "Y_LE", "Z_LE"]] = pd.DataFrame(
df["LEPoint"].tolist(), index=df.index
)
df[["X_TE", "Y_TE", "Z_TE"]] = pd.DataFrame(
df["TEPoint"].tolist(), index=df.index
)
fig = go.Figure()
fig.add_trace(
go.Scatter3d(
x=df["X_LE"],
y=df["Y_LE"],
z=df["Z_LE"],
mode="markers",
marker=dict(color="blue", size=8),
name="LEPoint",
)
)
fig.add_trace(
go.Scatter3d(
x=df["X_TE"],
y=df["Y_TE"],
z=df["Z_TE"],
mode="markers",
marker=dict(color="red", size=8),
name="TEPoint",
)
)
# Draw lines connecting LEPoint to TEPoint
for i in range(len(df)):
fig.add_trace(
go.Scatter3d(
x=[df["X_LE"][i], df["X_TE"][i]],
y=[df["Y_LE"][i], df["Y_TE"][i]],
z=[df["Z_LE"][i], df["Z_TE"][i]],
mode="lines",
line=dict(color="black", width=2),
showlegend=False,
)
)
# Setting same scale for all three axes
fig.update_layout(
scene=dict(aspectmode="cube", aspectratio=dict(x=1, y=1, z=1))
)
fig.show()
def draw_full(self, plot_lines=True):
plotly_fig = self.draw(backend="plotly")
# add line
if plot_lines:
df_line = self.data.lines_df
for i, row in df_line.iterrows():
plotly_fig.add_trace(
go.Scatter3d(
x=[row["@Point0"][0], row["@Point1"][0]],
y=[row["@Point0"][1], row["@Point1"][1]],
z=[row["@Point0"][2], row["@Point1"][2]],
mode="lines",
name="lines",
)
)
# ref line_application pt
plotly_fig.add_trace(
go.Scatter3d(
x=[self.data.lines_line_application_point[0]],
y=[self.data.lines_line_application_point[1]],
z=[self.data.lines_line_application_point[2]],
mode="markers",
marker=dict(
size=6,
),
)
)
# add pilot
plotly_fig.add_trace(
go.Scatter3d(
x=[self.data.pilot_pt[0]],
y=[self.data.pilot_pt[1]],
z=[self.data.pilot_pt[2]],
mode="markers",
marker=dict(
size=6,
),
)
)
plotly_fig.show()
def get_twist(self):
return [elt.twist for elt in self.wings[0].xsecs]
def _store_twist_origin(self):
self.data.origin_twist = self.get_twist()
def reset_twist(self):
self.modify_twist_airplane(self.data.origin_twist)
def modify_twist_airplane(self, twists: list):
"""
modify the twist without moving the LE point , just the angle
Args:
twists (list): The list of new twists in degger , abosolute to x axis.
Returns:
modify in place
"""
assert len(self.wings[0].xsecs) == len(
twists
), f"twists {len(twists)} isn't the same size as the number of sections of the airplane {len(self.wings[0].xsecs)}"
for i in range(len(twists)):
self.wings[0].xsecs[i].twist = twists[i]
def _compute_quasi_vertical_vect(self):
pt_center_aero = self.wings[0].aerodynamic_center()
pt_B = self.data.general["RisersDict"]["B"]
pt_B[1] = 0
return vector_between2pts(pt_center_aero, pt_B, unitize=True)
def _compute_main_pts(self):
"""
Compute key points for the glider: Carabiner point, the Pilot point, and Center of Gravity (CG) point.
The following input are coming from Ozcad source file:
- RiserLength [m]
- PilotVerticalCentreOfGravity [m]
- HarnessDragArea [m2]
- PilotLoad [kg]
The following input are hard coded: glider_mass ( from ParagliderData.__init__ )
Note:
- pt_B : center point between the Top B risers pt ( with y=0)
- quasi_vertical_vect : vector from the aerodynamic_center of the wing to pt_B
- Carabiner Point (carabiner_pt): The middle point of the main harness carabiner. extending quasi_vertical_vect from pt_B by the risers length
- Pilot Point (pilot_pt): the pilot's cg. Extending quasi_vertical_vect from carabiner_pt by the PilotVerticalCentreOfGravity
- Center of Gravity (cg): the combined center of gravity of the glider and the pilot.
"""
quasi_vertical_vect = self._compute_quasi_vertical_vect()
# carabiner pt
translate_vect_top_bottom_riser = (
quasi_vertical_vect * self.data.general["RiserLength"]
)
pt_B = self.data.general["RisersDict"]["B"]
pt_B[1] = 0
self.data.carabiner_pt = translate_point(pt_B, translate_vect_top_bottom_riser)
# Pilot pt
translate_vect_carabiner_pilot = (
quasi_vertical_vect * self.data.general["PilotVerticalCentreOfGravity"]
)
pilot_pt = translate_point(
self.data.carabiner_pt, translate_vect_carabiner_pilot
)
self.data.pilot_pt = pilot_pt
self.data.pilot_scx = self.data.general["HarnessDragArea"]
# cg pt
pt_center_aero = self.wings[0].aerodynamic_center()
pilot_mass_kg = self.data.general["PilotLoad"]
cg_pt = weighted_barycentric_average(
[pt_center_aero, self.data.pilot_pt], [self.data.glider_mass, pilot_mass_kg]
)
self.data.cg = cg_pt
def compute_triming(self):
cg_pt = np.array(self.data.cg)
LE_pt = self.wings[0].xsecs[0].xyz_le
twist0 = -self.wings[0].xsecs[0].twist
chord0 = self.wings[0].xsecs[0].chord
TE_pt = find_point0_TE(LE_pt, twist0, chord0)
Chord_line = Line.from_points(LE_pt, TE_pt)
projected_point = Chord_line.project_point(cg_pt)
ratio = distance_between_2pts(LE_pt, projected_point) / chord0
if projected_point[0] < LE_pt[0]:
ratio = -ratio
return ratio
def main_ChordAngle(self):
return self.data.ribs_df["ChordAngle"][0]
def derotate_all(self):
if not self.data.general["hasbeenDerotated"]:
self._derotate_wing()
self._derotate_lines()
self._compute_main_pts() # TODO group in one derotate
self.data.general["hasbeenDerotated"] = True
else:
print("Impossible : Glider already derotated")
def _derotate_lines(self):
# derotate the appplication point
angle_derotate = self.data.general["original_angle"]
print(f"The aplication lines point will be derotate by {angle_derotate} deg ")
self.data.lines_line_application_point = rotate_pt(
self.data.lines_line_application_point, angle_derotate
)
# derotate_each lines
self.data.lines_df["@Point0"] = self.data.lines_df["@Point0"].apply(
rotate_pt, theta_deg=angle_derotate
)
self.data.lines_df["@Point1"] = self.data.lines_df["@Point1"].apply(
rotate_pt, theta_deg=angle_derotate
)
# derotate riser webs
for riser, riser_pt in self.data.general["RisersDict"].items():
self.data.general["RisersDict"][riser] = rotate_pt(riser_pt, angle_derotate)
def _derotate_wing(self, angle_target=0):
"""Rotate the glider around y axis to bring the main chord to colinear with vector x.
Positive , is pitch down
Args:
angle_derotate (flaot): angle in degrees
"""
assert angle_target == 0, "angle_target != 0 not implemented"
angle_derotate = angle_target - self.data.general["original_angle"]
print(f"The glider will be derotate by {angle_derotate} deg ")
self.data.ribs_df["LEPoint_derot"] = self.data.ribs_df["LEPoint"].apply(
rotate_pt, theta_deg=-angle_derotate
)
# df_yaxis = pd.DataFrame(
# [
# {"y_axis": self.wings[0]._compute_frame_of_WingXSec(i)[1]}
# for i in range(len(self.wings[0].xsecs))
# ]
# )
# self.data.ribs_df = pd.concat([self.data.ribs_df, df_yaxis], axis=1)
# self.data.ribs_df["y_angle_rad"] = self.data.ribs_df["y_axis"].apply(
# angle_between
# )
# self.data.ribs_df["y_angle_deg"] = np.degrees(self.data.ribs_df["y_angle_rad"])
# self.data.ribs_df["correction_twist"] = np.cos(
# self.data.ribs_df["y_angle_rad"]
# ) * np.radians(angle_derotate)
derotate_wing_xsecs = []
for i, row in self.data.ribs_df.iterrows():
xsec = asb.WingXSec(
xyz_le=row["LEPoint_derot"],
chord=row["Chord"],
twist=0,
airfoil=asb.Airfoil(
row["Name"], coordinates=extract_profile(row["Profile2D"])
),
)
derotate_wing_xsecs.append(xsec)
derotate_wing = asb.Wing(
name="Main",
xsecs=derotate_wing_xsecs,
symmetric=True,
is_main_wing=True,
is_horizontal_stabilizer=False,
is_vertical_stabilizer=False,
is_fuselage=False,
)
self.wings[0] = derotate_wing
# %% TEMP PLOT
# def draw_skeleton(self):
# df = self.data.ribs_df
# # Extracting X, Y, and Z coordinates for LEPoint
# df[["X_LE", "Y_LE", "Z_LE"]] = pd.DataFrame(df["LEPoint"].tolist(), index=df.index)
# # Extracting X, Y, and Z coordinates for TEPoint
# df[["X_TE", "Y_TE", "Z_TE"]] = pd.DataFrame(df["TEPoint"].tolist(), index=df.index)
# # Plotting with isometric scaling
# fig = go.Figure()
# # Scatter plot for LEPoint
# fig.add_trace(
# go.Scatter3d(
# x=df["X_LE"],
# y=df["Y_LE"],
# z=df["Z_LE"],
# mode="markers",
# marker=dict(color="blue", size=8),
# name="LEPoint",
# )
# )
# # Scatter plot for TEPoint
# fig.add_trace(
# go.Scatter3d(
# x=df["X_TE"],
# y=df["Y_TE"],
# z=df["Z_TE"],
# mode="markers",
# marker=dict(color="red", size=8),
# name="TEPoint",
# )
# )
# # Draw lines connecting LEPoint to TEPoint
# for i in range(len(df)):
# fig.add_trace(
# go.Scatter3d(
# x=[df["X_LE"][i], df["X_TE"][i]],
# y=[df["Y_LE"][i], df["Y_TE"][i]],
# z=[df["Z_LE"][i], df["Z_TE"][i]],
# mode="lines",
# line=dict(color="black", width=2),
# showlegend=False,
# )
# )
# # Setting same scale for all three axes
# fig.update_layout(scene=dict(aspectmode="cube", aspectratio=dict(x=1, y=1, z=1)))
# fig.show()
# %% IMPORT GLIDERS
if __name__ == "__main__":
para = Paraglider(glider_model_path, derotate=False)
# para.draw_full()
# %%
foldername = "Enzo3"
tab_to_load = os.path.join(os.getcwd(), "xflr5_file", foldername, "tabs_Enzo3.xlsx")
speed_range = 100 # mm
df_tab = pd.read_excel(tab_to_load)
df_tab.interpolate(limit_direction="both", inplace=True)
tab_to_consider = ["A_average", "B"]
df_tab["distance_tab_pct"] = df_tab[tab_to_consider].diff(axis=1)[
tab_to_consider[-1]
]
df_tab["distance_tab_mm"] = (df_tab["distance_tab_pct"] / 100) * df_tab["chord"]
df_tab["speed_travel"] = df_tab["system"] * speed_range
df_tab["angle_speed_move"] = np.degrees(
np.arctan2(df_tab["speed_travel"], df_tab["distance_tab_mm"])
)
if para.data.general["HasCentreCell"]:
first_row = df_tab.iloc[[0]].copy()
df_tab = pd.concat([first_row, df_tab], ignore_index=True)
df_tab["angle_tot"] = df_tab["angle_speed_move"] - para.data.origin_twist
ax_speed = df_tab["angle_speed_move"].plot()
ax_speed.set_title("angle_change (deg)")
# %
para.modify_twist_airplane(df_tab["angle_tot"])
# %%
vlm_TS = asb.VortexLatticeMethod(
chordwise_resolution=8,
spanwise_resolution=1,
airplane=para,
op_point=asb.OperatingPoint(
velocity=12, # m/s
alpha=3, # degree
),
)
vlm_enzo_TS_result = vlm_TS.run()
# %%
ctr_s = asb.ControlSurface(
"flapo", symmetric=True, deflection=10, hinge_point=0.76, trailing_edge=True
)
w = para.wings[0]
for i in range(10):
w.xsecs[i].control_surfaces.append(ctr_s)
para_d = para.with_control_deflections({"flapo": 20})
# %%
vlm_TS_d = asb.VortexLatticeMethod(
chordwise_resolution=8,
spanwise_resolution=1,
airplane=para_d,
op_point=asb.OperatingPoint(
velocity=12, # m/s
alpha=3, # degree
),
)
vlm_enzo_TS_result_d = vlm_TS_d.run()
# %%
# cg_pt = np.array(para.data.cg)
# LE_pt = para.wings[0].xsecs[0].xyz_le
# twist0 = -para.wings[0].xsecs[0].twist
# chord0 = para.wings[0].xsecs[0].chord
# TE_pt = find_point0_TE(LE_pt, twist0, chord0)
# Chord_line = Line.from_points(LE_pt, TE_pt)
# projected_point = Chord_line.project_point(cg_pt)
# ratio = distance_between_2pts(LE_pt, projected_point) / chord0
# if projected_point[0] < LE_pt[0]:
# ratio = -ratio
# # %%
# angle_radians = -np.radians()
# line_direction = np.array([np.cos(angle_radians), 0, np.sin(angle_radians)])
# # Calculate the vector from the line's point to the given point
# point_to_line_point = cg_pt - np.array(LE_pt)
# # Calculate the projection of the point onto the line
# projection = (
# np.array(LE_pt) + np.dot(point_to_line_point, line_direction) * line_direction
# )
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment