% Retrieve or compute parameters for ducted-fan VTOL micro-UAV.
%
% Copyright (C) 2024, Naoki Sean Pross, ETH Zürich
% This work is distributed under a permissive license, see LICENSE.txt

function [params] = uav_params()

% Unit of measurements unless specified otherwise are
% Mass in kg
% Lenghts in m
% Time in s
% Frequencies in Hz
% Angular velocities in rad / s
% Uncertainty / measurement errors in percentage (between 0 and 1)

params = struct();

% ------------------------------------------------------------------------
% Physical constants

params.physics = struct(...
  'Gravity', 9.81, ...
  'AirDensity', 1.204 ... % kg / m^3
);

% ------------------------------------------------------------------------
% Mechanical measurements

% Inertia marix at center of mass,
% numbers from CAD are in g * mm^2, convert to kg * m^2
J = [
  5.856E+08, 1.121E+06, -3.506E+05;
  1.121E+06, 4.222E+08, 6.643E+06;
  -3.506E+05, 6.643E+06, 4.678E+08
] * 1e-9;

% Approximate propeller with a disk
m_prop = 200e-3; % mass of propeller
r_prop = 85e-3; % radius of propeller
J_prop = .5 * m_prop * r_prop^2;

params.mechanical = struct(...
  'Mass', 3.5, ...
  'DuctRadius', 46e-3, ...
  'DuctHeight', 171e-3, ...
  'FlapZDistance', 98e-3, ... % flap distance along z from center of mass
  'InertiaTensor', J, ...
  'GyroscopicInertiaZ', J_prop, ... % assume small angle
  'GyroscopicInertiaZUncertainty', .01 ... % in %
);

% ------------------------------------------------------------------------
% Actuator limits and measurements

% Servos usually give a "speed" in seconds / 60 degrees without load
servo_speed = 0.13; % seconds / 60 deg
servo_angular_velocity = (60 / servo_speed) * (pi / 180); % rad / s

params.actuators = struct(...
  'PropellerMaxAngularVelocity', 620.7, ... % in rad / s
  'ServoAbsMaxAngle', 20 * pi / 180, ... % in radians
  'ServoMaxTorque', 4.3 * 1e-2, ... % in kg / m
  'ServoNominalAngularVelocity', servo_angular_velocity ...
);

% IMU runs in NDOF mode
params.measurements = struct(...
  'SensorFusionDelay', 30e-3, ... % in s
  'LIDARAccuracy', 10e-3, ... % in m
  'LIDARMaxDistance', 40, ... % inm
  'LIDARBandwidth', 100, ... % in Hz for z position
  'IMUBandwidth', 100 ... % in Hz
);

% ------------------------------------------------------------------------
% Aerodynamics modelling

% Compute thrust proportionality factor from actuator limits
% from thrust relation F = k * omega^2.
% FIXME: this is not ideal, need better measurements
omega_max = params.actuators.PropellerMaxAngularVelocity;
F_max = 38.637; % in N (measured)
k_T = F_max / omega_max^2;

% Lift coefficient from data
% flap_stair = readtable('meas/Flap_Force_Stair.csv');
% [~, istart] = min(abs(flap.Force));
% [~, iend] = min(abs(flap.Force - 3));

% FIXME: LiftCoefficient comes from
% https://scienceworld.wolfram.com/physics/LiftCoefficient.html
params.aerodynamics = struct(...
  'ThrustOmegaProp', k_T, ... % in s^2 / (rad * N)
  'ThrustOmegaPropUncertainty', .01, ... % in %
  'FlapArea', 23e-3 * 10e-3, ... % in m^2
  'FlapAreaUncertainty', .02, ... % in %
  'DragCoefficients', [0, 0], ... % TODO
  'DragCoefficientsUncertainties', [.1, .1], ... % in %
  'LiftCoefficient', 2 * pi, ... % TODO
  'LiftCoefficientUncertainty', .05 ...% in %
);


% ------------------------------------------------------------------------
% Linearization point of non-linear dynamics.

% Compute theoretical thrust required to make the UAV hover
% from the relation mg = k * omega^2
% FIXME: This value should probably be replaced with a measurement
g = params.physics.Gravity;
m = params.mechanical.Mass;
k = params.aerodynamics.ThrustOmegaProp;

omega_hover = sqrt(m * g / k);

params.linearization = struct(...
  'PadeApproxOrder', 3, ...
  'Position', [0; 0; -2], ... % in inertial frame, z points down
  'Velocity', [0; 0; 0], ... % in inertial frame
  'Angles', [0; 0; pi / 4], ...   % in body frame
  'AngularVelocities', [0; 0; 0], ... % in body frame
  'Inputs', [0; 0; 0; 0; omega_hover] ... % Flaps at rest and turbine at X
);

% ------------------------------------------------------------------------
% Normalization (maximum values)

params.normalization = struct(...
  'HPosition', 1, ... % m
  'VPosition', 1, ...% m
  'HSpeed', 2, ... % m / s
  'VSpeed', 2, ... % m / s
  'PitchRollAngle', 10 * pi / 180, ... % rad
  'YawAngle', 5 * pi / 180, ... % rad
  'AngularRate', 2 * pi / 180, ... % rad / s
  'WindSpeed', 2 ... % m / s
);
end
% vim: ts=2 sw=2 et: