function y = fcn(xin,uin)    %12个状态 4输入 12输出  x1滚转角 x3俯仰角 x5偏航角 x7高度z x9平动x x11平动y
u=uin(1:4);
x=xin(1:12);
Dx=10^-6;
Dy=10^-6;
Dz=10^-4;
g=9.81;
m=1;
Jx=8.1*10^-3;
Jy=8.1*10^-3;
Jz=14.2*10^-3;
jr=1.04*10^-4;
l=0.24;
k01=5.42*10^-5;
k02=1.1*10^-6;
%参数
R=[k01 k01 k01 k01;-k01 k01 k01 -k01;-k01 -k01 k01 k01;-k02 k02 -k02 k02];
omega_2=R^-1*u;
omega1=sqrt(omega_2(1));
omega2=sqrt(omega_2(2));
omega3=sqrt(omega_2(3));
omega4=sqrt(omega_2(4));   %转速
uf=-omega1+omega2-omega3+omega4; 
dx1=x(2);                                                      %滚转角一阶导
dx2=((Jy-Jz)/Jx)*x(4)*x(6)+(jr/Jx)*x(4)*uf+(l/Jx)*u(2);     %滚转角二阶导
dx3=x(4);                                                      %俯仰角一阶导
dx4=((Jz-Jx)/Jy)*x(2)*x(6)-(jr/Jy)*x(2)*uf+(l/Jy)*u(3);     %俯仰角二阶导
dx5=x(6);                                                      %偏航角一阶导
dx6=((Jx-Jy)/Jz)*x(2)*x(4)+(1/Jz)*u(4);                     %偏航角二阶导
dx7=x(8);                                                      %高度z一阶导
dx8=(1/m)*u(1)*cos(x(3))*cos(x(1))-g-(1/m)*Dz*x(8)*x(8);    %高度z二阶导
dx9=x(10);                                                      %x一阶导
dx10=(1/m)*u(1)*(sin(x(3))*cos(x(5))*cos(x(1))+sin(x(5))*sin(x(1)))-(1/m)*Dx*x(10)*x(10);%x二阶导
dx11=x(12);                                                      %y一阶导
dx12=(1/m)*u(1)*(sin(x(3))*sin(x(5))*cos(x(1))-sin(x(1))*cos(x(5)))-(1/m)*Dy*x(12)*x(12);%y二阶导
y =[dx1;dx2;dx3;dx4;dx5;dx6;dx7;dx8;dx9;dx10;dx11;dx12];