function y = obs_TC04c_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)

% Constants
d2r   = pi/180;
r2d   = 180/pi;

% Geometry data
G0     =    9.80665D+0;

% 運動方程式まわりの設定(運動方程式をどの軸で設定したかによる)
global stables;
U0     = stables(1);  % トリム状態における対気速度[m/s]
W0     = stables(2);  % トリム状態における下方向速度[m/s]
THETA0 = stables(3);  % トリム状態におけるピッチ角

% State variables
U     = x(1);
Alpha = x(2);
Theta = x(3);
Q     = x(4);

% Input variables
de    = u(1);
dt    = zeros(size(u(1)));

% Parameters
Xu   = param( 1);
Xa   = param( 2);
Xdt  = 0;
Zu   = param( 3);
Za   = param( 4);
Zq   = param( 5);
Zde  = param( 6);
Zdt  = 0;
Mu   = param( 7);
Ma   = param( 8);
Mq   = param( 9);
Mde  = param(10);
Mdt  = 0;

Alpha0 = W0 / U0;
W = Alpha * U0;

% 運動方程式より微分量
Udot = Xu * U + Xa * Alpha - W0 * Q - G0 * cos(THETA0) * Theta + Xdt * dt;

%Aldot = (Zu * U + Za * Alpha + (U0 + Zq) * Q - G0 * sin(THETA0) * Theta ...
%      + Zde * de + Zdt * dt) / U0;

%Wdot = Aldot * U0;

Wdot = Zu * U + Za * Alpha + (U0 + Zq) * Q - G0 * sin(THETA0) * Theta ...
      + Zde * de + Zdt * dt;

% Observations
% Z 04 TASCG[M/S]   真大気速度(C.G.)  => U/W_{0}, u/w から生成
% Z 05 ALFCG[RAD]   迎角(C.G.)        => alpha
% Z 06 THE[RAD]     ピッチ(Theta)     => theta0 + theta
% Z 07 Q[RAD/S]     ピッチ角速度      => q
%   08 QDOT[RAD/S2] ピッチ角加速度    => qから運動方程式で生成
% Z 09 AXCG[M/S2]   加速度X(C.G.)     => uから運動方程式で生成
% Z 10 AZCG[M/S2]   加速度Z(C.G.)     => dot{w} = U_{0} * dot{alpha}
y(1)  =  ((U + U0)^2 + (W + W0)^2)^0.5;
y(2)  =  Alpha0 + Alpha;
y(3)  =  THETA0 + Theta;
y(4)  =  Q;

cos_y3 = cos(y(3));
sin_y3 = sin(y(3));

y(5)  = -Udot - Q * (W0 + W) + G0 * sin_y3; 
y(6)  =  Wdot - Q * (U0 + U) - G0 * cos_y3;

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

return
% end of function