"""
The ``Trajectory`` is the interface to the differential equation solver
for the disc trajectory.
"""
from dataclasses import dataclass
from numbers import Number
from typing import Dict, Union
import numpy as np
[docs]def rotation_matrix(phi: float, theta: float) -> np.ndarray:
"""
Compute the (partial) rotation matrix that transforms from the
lab frame to the disc frame. Note that because of azimuthal
symmetry, the azimuthal angle (`gamma`) is not used.
"""
sp, cp, st, ct = (np.sin(phi), np.cos(phi), np.sin(theta), np.cos(theta))
return np.array(
[[ct, sp * st, -st * cp], [0, cp, sp], [st, -sp * ct, cp * ct]]
)
[docs]@dataclass
class Trajectory:
"""
Class for computing the disc flight trajectory. Takes initial values
and interfaces with an ODE solver.
Units are meters [m] for length, kilograms [kg] for mass, seconds [s]
for time, and radians [rad] for angles.
Args:
x (float): horizontal position; default is 0 m
y (float): horizontal position; default is 0 m
z (float): vertical position; default is 1 m
vx (float): x-velocity; default is 10 m/s
vy (float): y-velocity; default is 0 m/s
vz (float): z-velocity; default is 0 m/s
phi (float): 1st Euler angle (pitch); default is 0 rad
theta (float): 2nd Euler angle (roll); default is 0 rad
gamma (float): 3rd Euler angle (spin); default is 0 rad
phidot (float): phi angular velocity; default is 0 rad/s
thetadot (float): theta angular velocity; default is 0 rad/s
gammadot (float): gamma angular velocity; default is 50 rad/s
"""
x: float = 0
y: float = 0
z: float = 1
vx: float = 10
vy: float = 0
vz: float = 0
phi: float = 0
theta: float = 0
gamma: float = 0
phidot: float = 0
thetadot: float = 0
gammadot: float = 50
[docs] def __post_init__(self):
self._initial_conditions: Dict[str, float] = {
coord: getattr(self, coord)
for coord in [
"x",
"y",
"z",
"vx",
"vy",
"vz",
"phi",
"theta",
"gamma",
"phidot",
"thetadot",
"gammadot",
]
}
[docs] def reset(self) -> None:
for k, v in self._initial_conditions.items():
setattr(self, k, v)
return
@property
def velocity(self) -> np.ndarray:
return np.array([self.vx, self.vy, self.vz])
@velocity.setter
def velocity(self, v: np.ndarray) -> None:
self.vx, self.vy, self.vz = v
return
@property
def angular_velocity(self) -> np.ndarray:
return np.array([self.phidot, self.thetadot, self.gammadot])
@angular_velocity.setter
def angular_velocity(self, v: np.ndarray) -> None:
self.phidot, self.thetadot, self.gammadot = v
return
[docs] def derived_quantities(
self,
) -> Dict[str, Union[float, np.ndarray, Dict[str, np.ndarray]]]:
"""
Compute intermediate quantities on the way to computing the time
derivatives of the kinematic variables.
"""
# Rotation matrix
R = rotation_matrix(self.phi, self.theta)
# Unit vectors
zhat = R[2]
velocity = self.velocity
v_dot_zhat = velocity @ zhat
v_in_plane = velocity - zhat * v_dot_zhat
xhat = v_in_plane / np.linalg.norm(v_in_plane)
yhat = np.cross(zhat, xhat)
# Angle of attack
angle_of_attack = -np.arctan(v_dot_zhat / np.linalg.norm(v_in_plane))
# Disc angular velocities
angular_velocity = self.angular_velocity
w_prime = np.array(
[
angular_velocity[0] * np.cos(self.theta),
angular_velocity[1],
angular_velocity[0] * np.sin(self.theta) + angular_velocity[2],
]
)
# Angular velocity in lab coordinates
w_lab = w_prime @ R
# Angular velocity components along the unit vectors in the lab frame
U = np.array([xhat, yhat, zhat])
w = U @ w_lab
return {
"unit_vectors": {"xhat": xhat, "yhat": yhat, "zhat": zhat},
"angle_of_attack": angle_of_attack,
"rotation_matrix": R,
"w_prime": w_prime,
"w_lab": w_lab,
"w": w,
}