function [dxdt] = Function(x,Aa,Ab,V0A,V0B,Rg,Ts,F_LVDT,m_T,Pa,Ps,cd_ISO,V,b,bS,time,t,vs,c1,c2,Fc1,Fc2,delta,toll_v,Fs,s,l1,coeff)

global a

V_mod=interp1(time,V,t);

if V_mod>5 
    A_eqmanS=abs((-4.929e-07)*(10-V_mod)+2.372e-06); 
    cond=A_eqmanS<0;
    A_eqmanS(cond)=0;
    cond=x(1)/Pa<=bS; %subsonico
    gA(cond)=-cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*x(1)*A_eqmanS*sqrt(1-((x(1)/Pa-bS)/(1-bS))^2); 
    cond=x(1)/Pa>bS; %sonico
    gA(cond)=-cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*x(1)*A_eqmanS; 
    % 4 CARICO
    A_eqman=abs((3.869e-07)*V_mod-1.998e-06);
    cond=A_eqman<0;
    A_eqman(cond)=0;
    cond=x(2)/Ps>=b; %subsonico
    gB(cond)=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*A_eqman*sqrt(1-((x(2)/Ps-b)/(1-b))^2);
    cond=x(2)/Ps<b; %sonico
    gB(cond)=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*A_eqman;
    
elseif V_mod==5 
    gA=0;
    gB=0;

else 
    A_eqman=abs((-3.953e-07)*V_mod+1.913e-06);
    cond=A_eqman<0;
    A_eqman(cond)=0;
    cond=x(1)/Ps>=b; %subsonico
    gA(cond)=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*A_eqman*sqrt(1-((x(1)/Ps-b)/(1-b))^2);
    cond=x(1)/Ps<b; %sonico
    gA(cond)=cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*Ps*A_eqman;
    % 4 SCARICO
    A_eqmanS=abs((5.084e-07)*(10-V_mod)-2.637e-06);
    cond=A_eqmanS<0;
    A_eqmanS(cond)=0;
    cond=x(2)/Pa<=bS; %subsonico
    gB(cond)=-cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*x(2)*A_eqmanS*sqrt(1-((x(2)/Pa-bS)/(1-bS))^2);
    cond=x(2)/Pa>bS; %sonico
    gB(cond)=-cd_ISO*sqrt(293/Ts)*0.686/sqrt(Rg*Ts)*x(2)*A_eqmanS; 

end

GA=gA*coeff;
GB=gB*coeff;


Fp=(x(1)-Pa)*Aa-(x(2)-Pa)*Ab;
v=x(4);
%1)Caso di velocità nulla
if abs(v)-toll_v<=0  
cond=abs(Fp-F_LVDT)<Fs;
Ff(cond)=Fp-F_LVDT; 
cond=abs(Fp-F_LVDT)>=Fs;
Ff(cond)=Fs*sign(Fp-F_LVDT);

else
%2a) Velocità crescenti
cond=a>=0;
Ff(cond)=c1*v+(Fc1+(Fs-Fc1)*exp(-(abs(v)./vs).^delta)).*sign(v);
%2a) Velocità decrescenti
cond=a<0;  
Ff(cond)=c2*v+Fc2*sign(v);


end
F_att=(Ff);

dxdt(1,1) = ((Rg*Ts)/(Aa*(x(3)+l1)+V0A))*(GA-Aa*x(1)*x(4)/(Rg*Ts));
dxdt(2,1) = ((Rg*Ts)/(Ab*(s-x(3)-l1)+V0B))*(GB+Ab*x(2)*x(4)/(Rg*Ts));
dxdt(3,1) = x(4);
dxdt(4,1) = 1/m_T*(Fp-F_att-F_LVDT);

a=dxdt(4,1);
end