function out=FOC(udc,isa,isb,spdfeed)
global Rs Rr ts np Lm Lr Ls sigma;
global PWM;
global Fsalfa Fsbeta;
%######################################
%rebuild usalfa, usbeta by sa,sb,sc,udc with none-deadzone
%######################################
sa=PWM(1);
sb=PWM(2);
sc=PWM(3);
usalfa=sqrt(2/3)*udc*(sa-1/2*sb-1/2*sc);
usbeta=sqrt(2/3)*udc*(sqrt(3)/2*sb-sqrt(3)/2*sc);
%######################################
%stator current 3s-2s transform
%######################################
isc=-isa-isb;
isalfa=(isa-1/2*isb-1/2*isc)*2/3;
isbeta=(sqrt(3)/2*isb-sqrt(3)/2*isc)*2/3;
%######################################
%U-I model, caculate the rotor flux
%######################################
% EMFalfa=usalfa-isalfa*Rs;
% EMFbeta=usbeta-isbeta*Rs;
% % FsAMP=sqrt(Fsalfa*Fsalfa+Fsbeta*Fsbeta);
%
% Fsalfa=Fsalfa+EMFalfa*ts;
% Fsbeta=Fsbeta+EMFbeta*ts;
%
% Fralfa=Lr/Lm*(Fsalfa-sigma*Ls*isalfa);
% Frbeta=Lr/Lm*(Fsbeta-sigma*Ls*isbeta);
%
%
%
% if Fralfa==0&Frbeta>0
% Frtheta=pi/2;
% elseif Fralfa==0&Frbeta==0
% Frtheta=0;
% elseif Fralfa==0&Frbeta<0
% Frtheta=3*pi/2;
% elseif Fralfa>0&Frbeta>=0
% Frtheta=atan(Frbeta/Fralfa);
% elseif Fralfa>0&Frbeta<0
% Frtheta=atan(Frbeta/Fralfa)+2*pi;
% else
% Frtheta=atan(Frbeta/Fralfa)+pi;
% end
global imr Frtheta;
cosine=cos(Frtheta);
sine=sin(Frtheta);
%build flux first
% if FsAMP<0.5
% iq=-sine*isalfa+cosine*isbeta;
% ud_ref=0;
% uq_ref=iqCTRL(iq,10,220,5,0.001,0);
% ualfa=cosine*ud_ref-sine*uq_ref;
% ubeta=sine*ud_ref+cosine*uq_ref;
% PWM=svpwm(ualfa,ubeta);
% out=[ualfa ubeta];
% return;
% end
id=cosine*isalfa+sine*isbeta;
iq=-sine*isalfa+cosine*isbeta;
%theta角计算
imr=imr+ts*Rr/Lr*(id-imr); %稳态磁化电流
ws=Rr/Lr*iq/id;
wr=spdfeed*2*np*pi/60;
w1=wr+ws; %w1=wr+ws
Frtheta=Frtheta+w1*ts;
if Frtheta>2*pi
Frtheta=Frtheta-2*pi;
elseif Frtheta<-2*pi
Frtheta=Frtheta+2*pi;
end
iqref=spdCTRL(spdfeed,1000,25,1,0.01,0);
%M轴（d轴）给定电流有电机参数计算的到
ud_ref=idCTRL(id,5.25,220,5,0.1,0);
uq_ref=iqCTRL(iq,iqref,220,20,0.1,0);
%######################################
%2r/2s inv_park
%######################################
ualfa=cosine*ud_ref-sine*uq_ref;
ubeta=sine*ud_ref+cosine*uq_ref;
out=[ualfa ubeta Frtheta];
