function [sys,X0]=dcmotor(t,x,u,FLAG);

% DCMOTOR 
%       S-function describing the non-linear dynamic equations of a
%       simple dcmotor given by eqs. 3.3-1, 3.3-5 and 3.3-6 of the
%       textbook "Electromechanical Motion Devices" P. Krause, McGraw
%       Hill 1989.
%
%	SYS=dcmotor(t,x,u,FLAG) returns, depending on FLAG, certain
%	system values given the current values of time, t, state 
%       vector, x, and input vector, u.
%
%	FLAG is used to indicate the value to be returned in SYS.
%       For a continuous system the following flags are used:
%
%         FLAG=0 to initialize the system 
%	  FLAG=1 to return state derivatives
%         FLAG=3 to return system outputs 
%
% Calling DCMOTOR with a FLAG of zero: [SIZES]=DCMOTOR([],[],[],0),  
%       returns a vector, SIZES, which contains the size of the state 
%       vector and other parameters (see FLAG=0) as indicated below.
%         SIZES(1) number of states
%         SIZES(2) number of discrete states
%         SIZES(3) number of outputs
%         SIZES(4) number of inputs
%         SIZES(5) number of roots (currently unsupported)
%         SIZES(6) direct feedthrough flag
%         SIZES(7) number of sample times
%
% For more information and other options see SFUNC.

% Copyright JCCK, Jan 1999

% ------------------------------------------------------------------------ %

%% Here, you must declare all the parameters of the system as global
%  of course this means that they should be loaded in the MATLAB 
%  workspace before running this S-function

global Rf Ra Laa Lff Laf Bm Jm

if FLAG == 0,	% initialize system
  % SYS and X0
   % load system parameters
    dcm_par;
   % This system has 3 states (x,v,i), 1 input (e), and 3 output (x,v,i)
    nx = 3; % # of states
    ny = 3; % # of outputs
    nu = 3; % # of inputs
    sys = [nx; 0; ny; nu; 0; 0];
   % Initial Conditions
    X0 = zeros(nx,1);
elseif FLAG == 1, % states derivatives
  % 1. Here put State depended Parameters if any or else...
    If = x(1); % field current
    Ia = x(2); % armature current
    Wr = x(3); % rotor angular velocity
    Vf = u(1); % field winding voltage
    Va = u(2); % armature winding voltage
    Tl = u(3); % load torque disturbance
  % 2. Here write the State Equations.
    dx(1) = -Rf/Lff*If + 1/Lff*Vf;
    dx(2) = -Ra/Laa*Ia -Laf/Laa*If*Wr + 1/Laa*Va;
    dx(3) = -Bm/Jm*Wr + Laf/Jm*If*Ia - 1/Jm*Tl;
  sys = dx;
elseif FLAG == 3,  % system outputs
    y(1) = x(1);
    y(2) = x(2);
    y(3) = x(3);
  sys = y; 
else,
  sys = [];
end