%% TIME DOMAIN SOLVER-FUNCTION %%

function xdot = nl_ae_solver(t,x,V)
global Params
global i
global uc
global tc
global ctr
global case_no
% global beta

load Params_all

rho=1.23;
b=Params.b;
aa=Params.a;
s=Params.s;
mw=Params.Wing_Mass;
rcg=Params.Wing_cg;
moff=mw*rcg;
ka=Params.Pitch_stiffness;
kh=Params.Plunge_stiffness;
kal=Params.Pitch_stiffness_lin;
mh=Params.Plunge_mass;%-Params.Plunge_mass_e;
Ia=Params.Pitch_It-Params.Pitch_It_e;
pitch_dr=Params.Pitch_dr;
plunge_dr=Params.Plunge_dr;

plunge_friction=Params.Plunge_Fk;

pitch_cc=2*sqrt(kal*Ia);
pitch_c=pitch_dr*pitch_cc;

plunge_cc=2*sqrt(kh*mh);
plunge_c=plunge_dr*plunge_cc;


CLA=Params.CLA;
CMA=Params.CMA;
CLB=0;%Params.CLB;
CMB=0;%Params.CMB;

u=0;

%% Matrices elements 
h =x(1);
h_dot=x(2);
a  =x(3);
a_dot=x(4);

% Damping Loads
% FC_h=plunge_friction*tanh(10000*h_dot); %(activate for Coulomb Damping) 
FC_a=polyval([pitch_c 0],a_dot);
FC_h=polyval([plunge_c 0],h_dot);

% FC_h=0;
% FC_a=0;

% Stiffness Loads
FK_a=polyval(ka,a);
% FK_a=polyval([kal 0],a);
FK_h=polyval([kh 0],h);

% Aerodynamic Loads
if(V==0)
Mo=0; Lf=0;
else
Mo = -rho*b^2*(pi)*((0.5-aa)*V*b*a_dot)*s + rho*V^2*b^2*s*CMB*u ...
+(CMA*rho*V^2*b^2*s)*(a + (0.5 - aa)*b*a_dot/V+h_dot/V);
   
Lf   = rho*V^2*b*s*CLA*(a + h_dot/V + (0.5 - aa)*b*a_dot/V)  ...
    +  rho*b^2*pi*(V*a_dot)*s  +  rho*V^2*b*s*CLB*u;
end

M_a_a=0;%- pi*rho*b^4*(1/8-aa/2)*s;
M_a_h=0;%- (pi*rho*b^3)/2*s;
M  = [mh moff ;moff-M_a_h  Ia-M_a_a];
%% Representation of equations
x = x(:);
xdot([1 3],1) = x([2 4],1);
xdot([2 4],1) = -inv(M)*[FC_h; FC_a] - inv(M)*[FK_h; FK_a] + inv(M)*[-Lf;Mo];
% sprintf(['Run No=',num2str(ctr),'    case_no=',num2str(case_no),'    t=',num2str(t)])
