function xdot = xdot_TC04_hfb_lon(ts, x, u, param)

% Function to compute the state derivatives (i.e., RHS of system state equations) 
% test_case = 4 -- Longitudinal motion: HFB-320 Aircraft
% Nonlinear model in terms of non-dimensional derivatives as function of
% variables in the stability axes (V, alfa)   
%                  states  - V, alpha, theta, q
%                  outputs - V, alpha, theta, q, qdot, ax, az
%                  inputs  - de, thrust  
% State equations (5.86), (5.87)


% State variables
VT    = x(1);
Alfa  = x(2);
The   = x(3);
Qrate = x(4);

% Input variables
de    = u(1);
Fe    = u(2);

% Parameters
CD0  = param( 1);
CDV  = param( 2);
CDAL = param( 3);
CL0  = param( 4);
CLV  = param( 5);
CLAL = param( 6);
CM0  = param( 7);
CMV  = param( 8);
CMAL = param( 9);
CMQ  = param(10);
CMDE = param(11);

% Geometry data
G0     =    9.80665D+0;
S      =    2.4957D-1;  % 主翼面積 [m^2]
M      =    1.6000D+0;  % 重量[kg]
CBAR   =    222.48D-3;  % コード[m] @ M.A.C.
I_Y    =    0.1250D+0;  % I_Y 慣性モーメント[kg*m^2]
L_TX   =    0.1100D+0;  % 重心からスラストまでの距離(X) [m] c.g. = 665 mm from nose, thrust = 775 mm
L_TZ   =    0.0500D+0;  % 重心からスラストまでの距離(Z) [m]
V0     =    17.800D+0;  % トリム状態における対気速度[m/s]
SIGMAT =    0.0000D+0;  % スラスト軸の傾き[deg]
RHO    =    1.1250D+0;  % 大気密度[kg/m^3] @ 地表

% Intermediate variables
QBAR   = 0.5D0 * RHO *VT^2;                  

% Right sides of state equations (5.86)
VTdot = - S / M * QBAR * (CD0 + (CDV * VT / V0) + (CDAL * Alfa))...
    + G0 * sin(Alfa - The)...
    + Fe / M * cos(Alfa + SIGMAT);

Aldot = - S / M * QBAR / VT * (CL0 + (CLV * VT / V0) + (CLAL * Alfa))...
    + Qrate...
    + G0 / VT *cos(Alfa - The + SIGMAT)...
    - Fe / (M * VT) * sin(Alfa + SIGMAT);

Qdot = S * CBAR / I_Y * QBAR * (CM0 + (CMV * VT / V0) + (CMAL * Alfa) + (CMQ * Qrate * CBAR / V0 / 2) + CMDE * de) + Fe / I_Y * (L_TX * sin(SIGMAT) + L_TZ * cos(SIGMAT));

% State derivatives
xdot(1) = VTdot;
xdot(2) = Aldot;
xdot(3) = Qrate;
xdot(4) = Qdot;

% xdot must be a column vector
xdot = xdot';

return
% end of function