function  fm_exc(flag)
% FM_EXC define Automatic Voltage Regulators
%
% FM_EXC(FLAG)
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%       FLAG = 5 non-windup limits
%
%see also FM_EXCIN
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Version:   1.0.0
%
%E-mail:    fmilano@thunderbox.uwaterloo.ca
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

global Bus Exc Syn DAE Oxl Pss Cluster

type = Exc.con(:,2);
ty1 = find(type == 1);
ty2 = find(type == 2);
ty3 = find(type == 3);

vg = DAE.V(Exc.bus);
vrmax = Exc.con(:,3);
vrmin = Exc.con(:,4);
Td = Exc.con(:,10);
Tr = Exc.con(:,11);
A = Exc.con(:,12);
B = Exc.con(:,13);

if ty1
  vm_1  = DAE.x(Exc.vm(ty1));
  vr1_1 = DAE.x(Exc.vr1(ty1));
  vr2_1 = DAE.x(Exc.vr2(ty1));
  efd_1 = DAE.x(Exc.vf(ty1));
  m0 = Exc.con(ty1,5);
  T1 = Exc.con(ty1,6);
  T2 = Exc.con(ty1,7);
  T3 = Exc.con(ty1,8);
  T4 = Exc.con(ty1,9);
  K1 = m0.*T2./T1;
  K2 = m0 - K1;
  K3 = T3./T4;
  K4 = 1 - K3;
  vr = m0.*vr2_1 + K3.*(K1.*(Exc.vrif(ty1) - vm_1) + vr1_1);
end
if ty2
  vm_2  = DAE.x(Exc.vm(ty2));
  vr1_2 = DAE.x(Exc.vr1(ty2));
  vr2_2 = DAE.x(Exc.vr2(ty2));
  efd_2 = DAE.x(Exc.vf(ty2));
  Ka = Exc.con(ty2,5);
  Ta = Exc.con(ty2,6);
  Kf = Exc.con(ty2,7);
  Tf = Exc.con(ty2,8);
  K5 = Kf./Tf;
end
if ty3
  vm_3  = DAE.x(Exc.vm(ty3));
  vr3_3 = DAE.x(Exc.vr3(ty3));
  efd_3 = DAE.x(Exc.vf(ty3));
  Kr  = Exc.con(ty3,5);
  T2r = Exc.con(ty3,6);
  T1r = Exc.con(ty3,7);
  Kr1 = Kr.*T1r./T2r;
  Kr2 = Kr - Kr1;
  vf0 = Exc.con(ty3,8);
  v0  = Exc.con(ty3,9);
end

switch flag
 case 3

  % Updating of reference voltages
  Exc.vrif = Exc.vrif0;
  if Oxl.n, Exc.vrif(Oxl.exc) = Exc.vrif(Oxl.exc) - DAE.x(Oxl.v); end
  if Pss.n, Exc.vrif(Pss.exc) = Exc.vrif(Pss.exc) + DAE.x(Pss.vss); end
  if Exc.cluster
    Exc.vrif(Exc.cluster) = Exc.vrif(Exc.cluster) + DAE.x(Cluster.Vs(Cluster.exc));
  end

  DAE.f(Exc.vm)  = (vg - DAE.x(Exc.vm))./Tr;

  % differential equations
  if ty1
    DAE.f(Exc.vr1(ty1)) = (K2.*(Exc.vrif(ty1) - vm_1) - vr1_1)./T1;
    DAE.f(Exc.vr2(ty1)) = (K4.*(vr1_1 + K1.*(Exc.vrif(ty1) - vm_1)) ...
                           - m0.*vr2_1)./(T4.*m0);
    % windup limiter
    vr = min(vr,vrmax(ty1));
    vr = max(vr,vrmin(ty1));
    DAE.f(Exc.vf(ty1)) = ...
        (-efd_1+vr-ceiling(efd_1,A(ty1),B(ty1),1))./Td(ty1);
  end
  if ty2
    DAE.f(Exc.vr1(ty2)) = (Ka.*(Exc.vrif(ty2)-vm_2-vr2_2- ...
                                K5.*(efd_2))-vr1_2)./Ta;
    DAE.f(Exc.vr2(ty2)) = (-K5.*efd_2-vr2_2)./Tf;
    % non-windup limiter
    idx = find(vr1_2 > vrmax(ty2) & DAE.f(Exc.vr1(ty2)) > 0);
    DAE.f(Exc.vr1(ty2(idx))) = 0;
    idx = find(vr1_2 < vrmin(ty2) & DAE.f(Exc.vr1(ty2)) < 0);
    DAE.f(Exc.vr1(ty2(idx))) = 0;
    vr1_2 = min(vr1_2,vrmax(ty2));
    vr1_2 = max(vr1_2,vrmin(ty2));
    DAE.x(Exc.vr1(ty2)) = vr1_2;
    DAE.f(Exc.vf(ty2)) = (-efd_2 + vr1_2 - ...
                           ceiling(efd_2,A(ty2),B(ty2),1))./Td(ty2);
  end
  if ty3
    DAE.f(Exc.vr3(ty3)) = (Kr2.*(Exc.vrif(ty3) - vm_3) - vr3_3)./T2r;
    DAE.f(Exc.vf(ty3)) = 1000*((vr3_3+Kr1.*(Exc.vrif(ty3)-vm_3)+vf0).*(vg(ty3)./v0)-efd_3);
    % non-windup limiter
    idx = find(efd_3 >= vrmax(ty3) & DAE.f(Exc.vf(ty3)) > 0);
    if idx, DAE.f(Exc.vf(ty3)) = 0; end
    DAE.x(Exc.vf(ty3)) = min(DAE.x(Exc.vf(ty3)),vrmax(ty3));
    idx = find(efd_3 <= vrmin(ty3) & DAE.f(Exc.vf(ty3)) < 0);
    if idx, DAE.f(Exc.vf(ty3)) = 0; end
    DAE.x(Exc.vf(ty3)) = max(DAE.x(Exc.vf(ty3)),vrmin(ty3));
  end

 case 4

  % common AVR Jacobians
  DAE.Fx = DAE.Fx - sparse(Exc.vm,Exc.vm,1./Tr,DAE.n,DAE.n);
  DAE.Fy = DAE.Fy + sparse(Exc.vm,Bus.n+Exc.bus,1./Tr,DAE.n,2*Bus.n);

  if ty1
    % AVR jacobians
    DAE.Fx = DAE.Fx + sparse(Exc.vr1(ty1),Exc.vm(ty1),-K2./T1,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Exc.vr1(ty1),Exc.vr1(ty1),-1./T1,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Exc.vr2(ty1),Exc.vm(ty1),-K4.*K1./T4./m0,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Exc.vr2(ty1),Exc.vr1(ty1),K4./T4./m0,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Exc.vr2(ty1),Exc.vr2(ty1),-1./T4,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Exc.vf(ty1),Exc.vf(ty1), ...
                             -(1+ceiling(efd_1,A(ty1),B(ty1),2))./ ...
                             Td(ty1),DAE.n,DAE.n);
    idx = find(vr < vrmax(ty1) & vr > vrmin(ty1));
    if idx,
      DAE.Fx = DAE.Fx + sparse(Exc.vf(ty1(idx)),Exc.vr1(ty1(idx)), ...
                               K3(idx)./Td(ty1(idx)),DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Exc.vf(ty1(idx)),Exc.vm(ty1(idx)), ...
                               -K3(idx).*K1(idx)./Td(ty1(idx)),DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Exc.vf(ty1(idx)),Exc.vr2(ty1(idx)), ...
                               m0(idx)./Td(ty1(idx)),DAE.n,DAE.n);
    end
    for i = 1:length(ty1), excjac(Exc.syn(ty1(i)),Exc.vf(ty1(i))); end
  end

  if ty2
    % AVR jacobians
    DAE.Fx = DAE.Fx - sparse(Exc.vr1(ty2),Exc.vm(ty2),Ka./Ta,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx - sparse(Exc.vr1(ty2),Exc.vr1(ty2),1./Ta,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx - sparse(Exc.vr1(ty2),Exc.vr2(ty2),Ka./Ta,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx - sparse(Exc.vr1(ty2),Exc.vf(ty2),K5.*Ka./Ta,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx - sparse(Exc.vr2(ty2),Exc.vr2(ty2),1./Tf,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx - sparse(Exc.vr2(ty2),Exc.vf(ty2),K5./Tf,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx - sparse(Exc.vf(ty2),Exc.vf(ty2), ...
                             (1+ceiling(efd_2,A(ty2),B(ty2),2))./ ...
                             Td(ty2),DAE.n,DAE.n);
    % windup limiter
    idx = find(vr1_2 < vrmax(ty2) & vr1_2 > vrmin(ty2));
    if idx
      DAE.Fx = DAE.Fx + sparse(Exc.vf(ty2(idx)), ...
                               Exc.vr1(ty2(idx)),1./Td(ty2(idx)), ...
                               DAE.n,DAE.n);
    end
    for i = 1:length(ty2), excjac(Exc.syn(ty2(i)),Exc.vf(ty2(i))); end
  end

  if ty3
    % AVR jacobians
    DAE.Fx = DAE.Fx + sparse(Exc.vr3(ty3),Exc.vr3(ty3),-1./T2r,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Exc.vr3(ty3),Exc.vm(ty3),-Kr2./T2r,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Exc.vf(ty3),Exc.vf(ty3),-1000,DAE.n,DAE.n);
    idx = find(efd_3 < vrmax(ty3) & efd_3 > vrmin(ty3));
    if ~isempty(idx)
      DAE.Fx = DAE.Fx + sparse(Exc.vf(ty3(idx)),Exc.vr3(ty3(idx)), ...
                               1000.*vg(ty3(idx))./v0(idx),DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Exc.vf(ty3(idx)),Exc.vm(ty3(idx)), ...
                               -1000*Kr1(idx).*vg(ty3(idx))./v0(idx),DAE.n,DAE.n);
      DAE.Fy = DAE.Fy + ...
               sparse(Exc.vf(ty3(idx)),Bus.n+Exc.bus(ty3(idx)),1000*(vr3_3(idx)+Kr1(idx).* ...
               (Exc.vrif(ty3(idx))-vm_3(idx))+vf0(idx))./v0(idx),DAE.n,2*Bus.n);
      for i = 1:length(idx), excjac(Exc.syn(ty3(idx(i))),Exc.vf(ty3(idx(i)))); end
    end
  end

 case 5

  if ty2, fm_windup(Exc.vr1(ty2),vrmax(ty2),vrmin(ty2)); end
  if ty3, fm_windup(Exc.vf(ty3),vrmax(ty3),vrmin(ty3)); end
  
end

%===================================================================
function output = ceiling(vf,A,B,flag)
Se = A.*(exp(B.*abs(vf))-1);
switch flag
 case 1, output = Se.*vf;
 case 2, output = Se+A.*B.*abs(vf).*exp(B.*abs(vf));
end

%===================================================================
function excjac(h,k)
global Syn Exc DAE
% updating Syncronous Machine Jacobians
Td10 = Syn.con(h,11);
ord = Syn.con(h,5);
switch ord
 case 3,
  DAE.Fx(Syn.e1q(h),k) = 1/Td10;
 case 4,
  DAE.Fx(Syn.e1q(h),k) = 1/Td10;
 case 5.1,
  DAE.Fx(Syn.e1q(h),k) = 1/Td10;
 case 5.3,
  xd   = Syn.con(h,8);
  x1d  = Syn.con(h,9);
  xL   = Syn.con(h,6);
  DAE.Fx(Syn.e1q(h),k) = (xd+xL)/(x1d+xL)/Td10;
 otherwise,
  Td20 = Syn.con(h,12);
  Taa  = Syn.con(h,24);
  DAE.Fx(Syn.e1q(h),k) = (1-Taa/Td10)/Td10;
  DAE.Fx(Syn.e2q(h),k) = Taa/Td20/Td10;
end