function y = obs_TC04_hfb_lon(ts, x, u, param)

% Function to compute the observation variables (i.e., RHS of observation 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  
% Observation equations (5.88), (5.89), (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 * x(1)^2;

CD    =  CD0 + CDV*VT/V0 + CDAL*Alfa;
CL    =  CL0 + CLV*VT/V0 + CLAL*Alfa;
SALFA =  sin(Alfa);
CALFA =  cos(Alfa);
CX    =  CL*SALFA - CD*CALFA;
CZ    = -CL*CALFA - CD*SALFA;

% Observations
y(1)  =  VT;
y(2)  =  Alfa;
y(3)  =  The;
y(4)  =  Qrate;
y(5)  =  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));
y(6)  =  QBAR * S / M * CX + Fe / M * cos(SIGMAT);
y(7)  =  QBAR * S / M * CZ - Fe / M * sin(SIGMAT);

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

return
% end of function