## 2 Models and Methods
### Reference Frame and Equations of Motion of Seagulls  
Only the longitudinal motion is considered in the present paper; however, the seagull has three freedom ranges in the longitudinal direction, as shown in the coordinate system of Figure 1: horizontal displacement \(x_{\text{b}}\), vertical displacement \(z_{\text{b}}\), and pitch angle \(\theta\). The Newton-Euler-based equations of a seagull's body are simplified to:  
\[ \begin{array}{l}\frac{\text{d}u(t)}{\text{d}t}=\frac{1}{m}X(t)-g\sin\theta- q(t)w(t)\\ \frac{\text{d}u(t)}{\text{d}t}=\frac{1}{m}Z(t)+g\cos\theta+q(t)u(t)\\ I_{yy}\frac{\text{d}q(t)}{\text{d}t}=M(t)\\ \frac{\text{d}\theta(t)}{\text{d}t}=q(t)\end{array}  \tag{1}\]  
Figure 1: The coordinate system for seagulls. where \(\theta\) represents the pitch angle between the horizontal plane and the \(x_{\text{b}}\)-axis. The \(u\) and \(w\) represent \(x_{\text{b}}\) and \(z_{\text{b}}\) components of velocity, respectively. \(q\) represents the \(y_{\text{b}}\) component of angular velocity. \(X\) and \(Z\) are the resultant forces from wing roots along \(x_{\text{b}}\) and \(z_{\text{b}}\) axes, respectively, and \(M\) represents the resultant moment from wing roots around the \(y_{\text{b}}\)- axes. \(I_{\text{yy}}\) represents the pitching moment of inertia about \(y_{\text{b}}\)-axis. Set \(u_{0}(t)\), \(w_{0}(t)\), \(q_{0}(t)\), and \(\theta_{0}(t)\) as the solution for trimmed forward-flight of the seagull. They satisfy the following equation:  
\[ \begin{array}{l}\frac{\text{d}u_{0}(t)}{\text{d}t}=\frac{1}{m}X_{0}(t)-g\sin \theta_{0}(t)-q_{0}(t)w_{0}(t)\\ \frac{\text{d}u_{0}(t)}{\text{d}t}=\frac{1}{m}Z_{0}(t)-g\cos\theta_{0}(t)-q_{0 }(t)u_{0}(t)\\ I_{\text{yy}}\frac{\text{d}\theta_{0}(t)}{\text{d}t}=M_{0}(t)\\ \frac{\text{d}\theta_{0}(t)}{\text{d}t}=q_{0}(t)\end{array}  \tag{2}\]  
We can set that:  
\[ \begin{array}{l}u(t)=u_{0}(t)+\delta u(t)\\ w(t)=w_{0}(t)+\delta w(t)\\ q(t)=q_{0}(t)+\delta q(t)\\ \theta(t)=\theta_{0}(t)+\delta\theta(t)\\ X(t)=X_{0}(t)+\delta X(t)\\ Z(t)=Z_{0}(t)+\delta Z(t)\\ M(t)=M_{0}(t)+\delta M(t)\end{array}  \tag{3}\]  
Substituting Equations (3) and (4) into Equation (1), and neglecting the second- and higher-order terms, give:  
\[ \begin{array}{l}\frac{\text{d}\delta u(t)}{\text{d}t}=\frac{1}{m}\delta X(t )-g\delta\theta(t)\cos\theta_{0}(t)\\ -q_{0}(t)\delta w_{0}(t)-w_{0}(t)\delta q(t)\\ \frac{\text{d}\delta w(t)}{\text{d}t}=\frac{1}{m}\delta Z(t)-g\delta\theta(t) \sin\theta_{0}(t)\\ +q_{0}(t)\delta u(t)+u_{0}(t)\delta q(t)\\ I_{yy}\frac{\text{d}\delta q(t)}{\text{d}t}=\delta M(t)\end{array}  \tag{5}\]  
The disturbance values of the resultant forces and moment can be expressed as an infinite series according to [PERSON]'s theorem [13], and neglecting the second- and higher-order terms results in:  
\[ \begin{array}{l}\delta X(t)=\frac{\partial X(t)}{\partial u}\delta u(t)+ \frac{\partial X(t)}{\partial w}\delta w(t)+\frac{\partial X(t)}{\partial q }\delta q(t)\\ \delta Z(t)=\frac{\partial Z(t)}{\partial u}\delta u(t)+\frac{\partial Z(t)}{ \partial w}\delta w(t)+\frac{\partial Z(t)}{\partial q}\delta q(t)\\ \delta M(t)=\frac{\partial M(t)}{\partial u}\delta u(t)+\frac{\partial M(t)}{ \partial w}\delta w(t)+\frac{\partial M(t)}{\partial q}\delta q(t)\end{array}  \tag{6}\]  
The partial derivatives with respect to \(u\) are, respectively, defined as:  
\[ \begin{array}{l}\frac{\partial X(t)}{\partial u}=\lim_{\Delta u\to 0}\frac{ \Delta X(t)}{\Delta u}\approx\frac{\Delta X(t)}{\Delta u}\\ \frac{\partial Z(t)}{\partial u}=\lim_{\Delta u\to 0}\frac{\Delta Z(t)}{ \Delta u}\approx\frac{\Delta Z(t)}{\Delta u}\\ \frac{\partial M(t)}{\partial u}=\lim_{\Delta u\to 0}\frac{\Delta M(t)}{ \Delta u}\approx\frac{\Delta M(t)}{\Delta u}\end{array}  \tag{7}\]  
Substituting Equations (6) and (7) into Equation (5) gives:  
\[\left[ \begin{array}{l}\delta\dot{u}\\ \delta\dot{w}\\ \delta\dot{q}\\ \delta\dot{\theta}\end{array} \right]=\mathbf{A}(t)\left[ \begin{array}{l} \delta u\\ \delta w\\ \delta q\\ \delta\theta\end{array}