function [wr_hat_se,h_out] = aci_se(i_sq,i_sd,psi_rq,psi_rd,theta_psi_r,h_in,p)
% This function estimates the rotor speed of Induction machine
% by means of open-loop mathematic model.
%
% Inputs:  
%       i_sq = Stationary q-axis stator current (amp)
%       i_sd = Stationary d-axis stator current (amp)
%       psi_rq = Stationary q-axis rotor flux linkage (volt.sec/rad)
%       psi_rd = Stationary d-axis rotor flux linkage (volt.sec/rad)
%       theta_psi_r = Synchronously rotating rotor flux angle (rad)
%       h_in = [theta_psi_r_prev; We_prev]
%       p = [T; Rr; Lr; Lm; fc; diff_max_limit; diff_min_limit; fb]
% Output:
%       wr_hat_se = Estimated rotor electrically angular velocity (rad/sec)
%       h_out = [theta_psi_r_curr; We_curr]
% where 
%       T = sampling period (sec)
%       Rr = Rotor resistance referred to stator (ohm)
%       Lr = Rotor self inductance referred to stator (H)
%       Lm = Magnetizing inductance (H)
%       fc = cut-off frequency for low-pass filter (Hz) 
%       diff_max_limit = Maximum limit for differentiation (%)
%       diff_min_limit = Minimum limit for differentiation (%)
%       theta_psi_r_prev,theta_psi_r_curr = Previous and current rotor flux angle (rad)
%       We_prev,We_curr = Previous and current estimated syn. angular velocity (rad/sec)

% Defining constants based on machine parameters
Tr = p(3)/p(2);     % Tr = Lr/Rr, Rotor time constant; depending on Rr

% Computing the low-pass filter time constant
tc = 1/(2*pi*p(5));    % Tc = 1/(2*pi*fc)  (sec)

Wb = 2*pi*p(8);
   
% Constants    
K1_se = 1/(Wb*Tr);
K2_se = 1/(p(8)*p(1)); 
K3_se = tc/(tc+p(1)); 
K4_se = p(1)/(tc+p(1));

% Name the variables
i_qs_se = i_sq;
i_ds_se = i_sd;
psi_qr_se = psi_rq;
psi_dr_se = psi_rd;
theta_r_se = theta_psi_r;
theta_r_old = h_in(1);
wr_psi_r = h_in(2);

%  Slip computation
psi_r_2 = psi_dr_se*psi_dr_se + psi_qr_se*psi_qr_se;
w_slip = K1_se*(psi_dr_se*i_qs_se - psi_qr_se*i_ds_se)/psi_r_2;
   
% Synchronous speed computation  
if ((theta_r_se < p(6))&(theta_r_se > p(7)))
 		w_syn = K2_se*(theta_r_se-theta_r_old);
else  w_syn = wr_psi_r;
end       

% low-pass filter 
wr_psi_r = K3_se*wr_psi_r + K4_se*w_syn;
    
theta_r_old = theta_r_se;
wr_hat_se = w_syn - w_slip;

if (wr_hat_se > 1)
    wr_hat_se = 1;
elseif (wr_hat_se < -1)
     wr_hat_se = -1;      
end
 
% Return the results
h_out = [theta_r_se; wr_psi_r];