function aeroestim_alfadot_teta1(bCz_alfa,bCz_delta,bCm_alfa,bCm_delta,u,t_step,teta_i)

% function to compute (identify) the aerodynamical coefficients from flight 
% test of a system with the following dynamical model:
% dx/dt = A x(t) + B u(t) with the observable y(t) = C x(t) + D u(t)
% where:
% A=[bCz_alfa const1 1;0 0 1;bCm_alfa 0 0], const1=9.81*sin(teta_0)/u_0;
% B=[bCz_delta;0;bCm_delta];
% C=[bCz_alfa const1 1;0 1 0];
% D=[bCz_delta;0];

% this is a 2-D linearized version of a rigid body flight dynamics  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_0=286; % the initial stream velocity at infinity (relative to the flight system)  

alfa_0=(9.81/u_0)*(bCm_delta)/(-bCz_alfa*bCm_delta+bCm_alfa*bCz_delta);
teta_0=alfa_0;
const1=9.81*sin(teta_0)/u_0;

A1=[-10 const1 1;0 0 1;-15 0 0];B1=[-5;0;-10];C1=[-10 const1 1;0 1 0];D1=[-5;0];
m1=idss(A1,B1,C1,D1);
m1.As=[NaN const1 1;0 0 1;NaN 0 0];
m1.Bs=[NaN;0;NaN];
m1.Cs=[NaN const1 1;0 1 0];
m1.Ds=[NaN ;0];
m1.Ts=0;
m1.As(1,1)=m1.Cs(1,1);
m1.Bs(1,1)=m1.Ds(1,1);
m2=init(m1);

w=aero_sim(bCz_alfa,bCz_delta,bCm_alfa,bCm_delta,u,t_step,const1,teta_0,teta_i);
alfadot=w(:,1);
teta=w(:,2);
estimData=iddata([alfadot teta],u,t_step);
model=pem(estimData,m2)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute the simulated data required for identification

function w=aero_sim(bCz_alfa,bCz_delta,bCm_alfa,bCm_delta,u,t_step,const1,teta_0,teta_i)

A=[bCz_alfa const1 1;0 0 1;bCm_alfa 0 0];
B=[bCz_delta;0;bCm_delta];
C=[bCz_alfa const1 1;0 1 0];
D=[bCz_delta;0];
m=idss(A,B,C,D);
m.Ts=0;

init=[teta_i teta_0-teta_i-cos(teta_0)/sin(teta_0) 0]';

e=.1*randn(length(u),2);% construct an output noise

Data_noisy=iddata([],[u e],t_step);
m.noisevar=[.01 0;0 .01];

simData_noisy=sim(Data_noisy,m,init);
alfadot_noisy=simData_noisy.outputdata(:,1);
teta_noisy_i=simData_noisy.outputdata(:,2);
teta_noisy=teta_noisy_i+ones(length(teta_noisy_i),1)*(cos(teta_0)/sin(teta_0)+teta_0);
w=[alfadot_noisy teta_noisy_i ];

subplot(2,1,1)
plot(alfadot_noisy)
title('noisy alfadot with output white noise varvariance=1% vs time(seconds)')
xlabel('time(sec/100)')
ylabel('alfadot(radians)')

subplot(2,1,2)
plot(teta_noisy)
title('noisy teta with output white noise varvariance=1% vs time(seconds)')
xlabel('time(sec/100)')
ylabel('teta(radians)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute the noiseless data and the noisy-to-noiseless variances 

Data0=iddata([],u,t_step);
simData0=sim(Data0,m);
alfadot0=simData0.outputdata(:,1);
teta0=simData0.outputdata(:,2);

var_alfadot=var(alfadot_noisy-alfadot0)
var_teta=var(teta_noisy-teta0)