% Set plane and get data
%   iplane = 1 is 747
%   iplane = 2 is F-104

iplane = 2;
airdata;

% Set initial condition

v0 = [0.1; 0; 0; 0];

% Calculate M, Atilde, and A matrices

M = zeros(4,4);
Atilde = zeros(4,4);

M(1,1) = m*u0;
M(2,2) = Ix;
M(2,3) = -Ixz;
M(3,2) = -Ixz;
M(3,3) = Iz;
M(4,4) = 1.0;

Atilde(1,1) = Yb;
Atilde(1,2) = Yp;
Atilde(1,3) = -m*u0 + Yr;
Atilde(1,4) = m*g;

Atilde(2,1) = Lb;
Atilde(2,2) = Lp;
Atilde(2,3) = Lr;

Atilde(3,1) = Nb;
Atilde(3,2) = Np;
Atilde(3,3) = Nr;

Atilde(4,2) = 1.0;

% Calculate A matrix and eigendecomposition

A = M\Atilde;
[R,D] = eig(A);
lambda = diag(D);

% Calculate analytic solution
alpha = R\v0;

t = linspace(0,100,501);
w = exp(lambda*t);
w(1,:) = alpha(1)*w(1,:);
w(2,:) = alpha(2)*w(2,:);
w(3,:) = alpha(3)*w(3,:);
w(4,:) = alpha(4)*w(4,:);

v = R*w;

% Plot results
subplot(221);
plot(t,v(1,:));
title(sprintf('sideslip (beta) vs t for %s',stype));
xlabel('t (sec)');
ylabel('beta (degrees)');

subplot(222);
plot(t,v(2,:));
title(sprintf('roll rate (p) vs t for %s',stype));
xlabel('t (sec)');
ylabel('p (degrees/sec)');

subplot(223);
plot(t,v(3,:));
title(sprintf('yaw rate (r) vs t for %s',stype));
xlabel('t (sec)');
ylabel('r (degrees/sec)');

subplot(224);
plot(t,v(4,:));
title(sprintf('roll angle (phi) vs t for %s',stype));
xlabel('t (sec)');
ylabel('phi (degrees)');