## 4 Simulation Setup
### Aerodynamic Disturbance Model
#### 4.1.1 Wind Speed Profile Generation  
The wind velocity vector, which has a horizontal and vertical component, can be modeled as the sum of two vectors [28]:  
\[\mathbf{w}(t)=\begin{cases}w_{h}(t)\\ w_{v}(t)\end{cases}=\mathbf{w}_{0}+\Delta\mathbf{w}(t) \tag{29}\]  
with \(\mathbf{w}_{0}\) representing the static wind vector (which represents the mean value of the overall wind vector) and \(\Delta\mathbf{w}(t)\) representing the fast time-varying (or gust) component. This last component features higher frequencies and smaller amplitudes in comparison with the static wind vector. Since our goal is to test the performance of the system under a dynamic load, we neglect the static component.  
The behavior of the wind gust component, \(\Delta\mathbf{w}(t)\), follows the Dryden model, which defines the horizontal (\(\Delta w_{h}(t)\)) and vertical (\(\Delta w_{v}(t)\)) components of wind as a sum of sinusoidal components:  
\[\Delta w_{j}(t)=\sum_{i=1}^{n}a_{ji}sin\big{(}\Omega_{j,i}t+\varphi_{j,i}\big{)} \tag{30}\]  
where the generic subscript \(j\) is adopted in lieu of the subscripts \(h\) or \(v\). In Equation (30), \(\Omega_{j,i}\) and \(\varphi_{j,i}\) are randomly selected frequencies and phase shifts, \(n\) is the number of sinusoids, and \(a_{j,i}\) are the amplitudes of the sinusoids. For realistic results, it is suggested to consider \(\Omega_{j,i}\) between 0.1 and 1.5 rad/s. For simplicity, we select equally spaced frequencies in this range for both horizontal and vertical components (\(\Omega_{j,i}=\Omega_{i}\)). The amplitudes, \(a_{j,i}\), are given by the following equation:  
\[a_{j,i}=\sqrt{\Delta\Omega\Phi_{j}(\Omega_{i})} \tag{31}\]  
where \(\Delta\Omega\) is the interval between the chosen frequencies, and \(\Phi(\Omega_{i})\) is a Power Spectral Density (PSD) function evaluated for \(\Omega=\Omega_{i}\). The PSD expressions for the horizontal and vertical wind components are different and are defined in the Dryden model, respectively, as follows:\[\Phi_{h}(\Omega_{i})=\sigma_{h}^{2}\frac{2L_{h}}{\pi}\frac{1}{1+\left(L_{h} \Omega\right)^{2}} \tag{32}\] \[\Phi_{v}(\Omega_{i})=\sigma_{v}^{2}\frac{L_{v}}{\pi}\frac{1+3(L_{v} \Omega)^{2}}{\left(1+\left(L_{v}\Omega\right)^{2}\right)^{2}}\]  
where \(\sigma_{h}\) and \(\sigma_{v}\) are the horizontal and vertical turbulence intensities, and \(L_{h}\) and \(L_{v}\) are the horizontal and vertical gust length scales. As suggested in [29], the vertical length scale and the turbulence intensity can be assumed to be \(L_{v}=|z|\) (where \(z\) is the altitude of the UAV) and \(\sigma_{v}=0.1\omega_{20}\) (\(\omega_{20}\) is the given wind speed in knots at a 20 ft altitude). In [30], it is suggested to use \(\omega_{20}=15\) knots for light turbulence. Finally, \(L_{h}\) and \(\sigma_{h}\) are found with the following equations, where the parameters are modified to allow us to have the altitude input, \(z\), expressed in SI units:  
\[\begin{array}{l}\frac{L_{h}}{L_{v}}=\frac{1}{\left(0.177+0.0027z\right)^{1. 2}}\\ \frac{\sigma_{h}}{\sigma_{v}}=\frac{1}{\left(0.177+0.0027z\right)^{0.4}}\end{array} \tag{33}\]  
The altitude at which the UAM operates is set to \(z=3\) m.