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

% Function to compute the state derivatives (i.e., RHS of system state equations) 
% Lateral motion
% Nonlinear model in terms of dimensional derivatives as function of
% variables in the stability axes (V, alfa)   
%                  states  - beta, p, r, phi, psi
%                  outputs - beta, p, r, phi, psi, ay
%                  inputs  - da, dr

% Geometry data
G0     =    9.80665D+0;
global stables;
U0     = stables(1);  % トリム状態における対気速度[m/s]
W0     = stables(2);  % トリム状態における下方向速度[m/s]
THETA0 = stables(3);  % トリム状態におけるピッチ角(安定軸)

% State Variables
Bet = x(1);
Prate = x(2);
Rrate = x(3);
Phi = x(4);
Psi = x(5);

% Input variables
da    = u(1);
dr    = u(2);

% Parameters
Yb   = param( 1);
Yp   = param( 2);
Yr   = param( 3);
Ydr  = param( 4);
Lb   = param( 5);
Lp   = param( 6);
Lr   = param( 7);
Lda  = param( 8);
Ldr  = param( 9);
Nb   = param(10);
Np   = param(11);
Nr   = param(12);
Nda  = param(13);
Ndr  = param(14);

% Right sides of state equations (5.86)
Betdot = (Yb * Bet + (W0 + Yp) * Prate + (Yr - U0) * Rrate ...
    + G0 * cos(THETA0) * Phi + Ydr * dr) / U0;
Pdot = Lb * Bet + Lp * Prate + Lr * Rrate + Lda * da + Ldr * dr;
Rdot = Nb * Bet + Np * Prate + Nr * Rrate + Nda * da + Ndr * dr;
Phidot = Prate + Rrate * tan(THETA0);
Psidot = Rrate * sec(THETA0);

% State derivatives
xdot(1) = Betdot;
xdot(2) = Pdot;
xdot(3) = Rdot;
xdot(4) = Phidot;
xdot(5) = Psidot;

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

return
% end of function