function [sys,x0,str,ts]=motor(t,x,u,flag,Incond,Mcoff)
% A simulation s-function of 5-order motor math model,the syntax is   
%  [Ia;Ib;Ic;wr;theta]=motor(ua,ub,uc,TL,Mcoff)
% Where
%   [ ia, ib, ic ] is vector of motor stator phase current 
%   wr: angular velocity
%   theta:  angle of rotation
%   uc,ub,uc are three phases input stator volatges of induction motor.
%   TL is load torque
%   Mcoff is vector of motor cofficients, the definition format is:    Mcoff=[Rs,Lls,L0,Lms,Ld,Lq,p,J,pusa]
%   where
%       Rs: stator  resistance           
%       Lls:  leakage inductance
%       L0+Lls: the DC component of the phase inductance 
%       Lms: Peak value of the AC  component of the phase inductance 
%       Ld: direct component of the phase inductance
%       Lq: quadrature component of the phase inductance
%       p: pole-pair number
%       J:  inertia
%       pusa: Peak value of the permanent magnet flux linkage of the phase winding 

%   Designer:   Automation Institute, HUST, PRC
%   Version:     1.0
%   Date:         26, June, 2003

if nargin~=6
   error('Error in input arguments');
end
% The following outlines the general structure of an S-function.
switch flag,
  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,
    [sys,x0,str,ts]=mdlInitializeSizes(t,x,u,Incond);
  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1,
    sys=mdlDerivatives(t,x,u,Mcoff);
  %%%%%%%%%%%
  % Outputs %
  %%%%%%%%%%%
  case 3,
    sys=mdlOutputs(t,x,u,Mcoff);
    
case {2,4,9}
    sys=[];
  %%%%%%%%%%%%%%%%%%%%
  % Unexpected flags %
  %%%%%%%%%%%%%%%%%%%%
  otherwise
    error(['Unhandled flag = ',num2str(flag)]);
end
% end motor
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
function [sys,x0,str,ts]=mdlInitializeSizes(t,x,u,Incond)
% call simsizes for a sizes structure, fill it in and convert it to a
% sizes array.
% Note that in this example, the values are hard coded.  This is not a
% recommended practice as the characteristics of the block are typically
% defined by the S-function parameters.
sizes = simsizes;
% States    X=[ ia; ib; ic; wr; theta]
sizes.NumContStates  = 5;       
sizes.NumDiscStates  = 0;
% Output   Y=X
sizes.NumOutputs     = 6;        
% Input      U=[ ua; ub; uc; TL ]
sizes.NumInputs      = 4;         
sizes.DirFeedthrough = 0;
% at least one sample time is needed
sizes.NumSampleTimes = 1;   

sys = simsizes(sizes);
% initialize the initial conditions
x0  = Incond';
% str is always an empty matrix
str = [];
% initialize the array of sample times
ts  = [0 0];
%end mdlInitializeSizes
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%  Update state variables
function sys=mdlDerivatives(t,x,u,Mcoff)
% Initialize Motor Parameters
Rs=Mcoff(1);    Lls=Mcoff(2);   L0=Mcoff(3);    Lms=Mcoff(4);  Ld=Mcoff(5);   
Lq=Mcoff(6);    p=Mcoff(7);     J=Mcoff(8);    pusa=Mcoff(9);   Jtype=Mcoff(10);
% middle variables definition
wr=x(4);  
%wr=2*pi*50;

theta=x(5);
if(x(5)>=2*pi) 
    x(5)=theta-2*pi;
end
if(x(5)<=-2*pi)
    x(5)=theta+2*pi;
end
x(5)=theta;
twotheta=2*theta;

a0=-L0/2+Lms*cos(twotheta);  a1=-L0/2+Lms*cos(twotheta-0.666667*pi);   a2=-L0/2+Lms*cos(twotheta+0.666667*pi);

b0=-2*wr*Lms*sin(twotheta);  b1=-2*wr*Lms*sin(twotheta-0.666667*pi);   b2=-2*wr*Lms*sin(twotheta+0.666667*pi);

c0=pusa*wr*cos(theta);  c1=pusa*wr*cos(theta-0.666667*pi);   c2=pusa*wr*cos(theta+0.666667*pi);

A=[ Lls+1.5*L0+a0,  a1, a2; a1, Lls+1.5*L0+a2,  a0; a2, a0, Lls+1.5*L0+a1   ];

B=[ Rs+b0,  b1, b2; b1, Rs+b2,  b0; b2, b0, Rs+b1   ];

C=[ c0; c1; c2  ];

if (cond(A)==inf)      
     %warning(' Matrix A is a singular matrix' ) ; 
     error(' Matrix A is a singular matrix, the simulation is halted' ) ; 
end 
 Ainv=inv(A);
 xnew=[ x(1); x(2); x(3)];
 unew=[ u(1); u(2); u(3)];
sysone = Ainv*(-B*xnew-C+unew);

E1=[ 1,  -0.5,   -0.5;   -0.5,   -0.5,   1;  -0.5,   1,  -0.5 ];
E2=[ 0,  -1,   1;   -1,   1,   0;  1,   0,  -1 ];
E3=[ 1,  -0.5,   -0.5 ];
E4=[ 0,  1,   -1 ];

Te=p*(Ld-Lq)*0.333333*( (xnew'*E1*xnew)*sin(twotheta) +0.866*(xnew'*E2*xnew)*cos(twotheta) )+...
    p*pusa*(  E3*xnew*cos(theta)+0.866*E4*xnew*sin(theta)  );

systwo=[ (Te-u(4))*p/J; wr];

sys = [sysone; systwo];
% end mdlDerivatives

%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%
function sys=mdlOutputs(t,x,u,Mcoff)
p=Mcoff(7);Ld=Mcoff(5);Lq=Mcoff(6); pusa=Mcoff(9); 
E1=[ 1,  -0.5,   -0.5;   -0.5,   -0.5,   1;  -0.5,   1,  -0.5 ];
E2=[ 0,  -1,   1;   -1,   1,   0;  1,   0,  -1 ];
E3=[ 1,  -0.5,   -0.5 ];
E4=[ 0,  1,   -1 ];

xnew=[ x(1); x(2); x(3)];
theta=x(5);
if(x(5)>=2*pi) 
    x(5)=theta-2*pi;
end
if(x(5)<=-2*pi)
    x(5)=theta+2*pi;
end
x(5)=theta;
twotheta=2*theta;

Te=p*(Ld-Lq)*0.333333*( (xnew'*E1*xnew)*sin(twotheta) +0.866*(xnew'*E2*xnew)*cos(twotheta) )+...
    p*pusa*(  E3*xnew*cos(theta)+0.866*E4*xnew*sin(theta)  );

x4new=9.5493*x(4)/p;
sys = [ xnew;x4new;x(5);Te ];

% end mdlOutputs