function fm_sae1(flag)% FM_SAE1 define a subtransmission area equivalent%         with three loads and three LTCs. This model%         uses one equivalent state variables presenting%         a "mean" value of the LTCs' tap ratios.%% FM_SAE1(FLAG)%       FLAG = 0 initialization%       FLAG = 1 algebraic equations%       FLAG = 2 algebraic Jacobians%       FLAG = 3 differential equations%       FLAG = 4 state matrix%       FLAG = 5 non-windup limits%%see also FM_SAE2 and FM_SAE3%%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 Milanoglobal Bus DAE SAE1for i = 1:SAE1.n   if isempty(SAE1.m)      m = 1;   else      m = DAE.x(SAE1.m(i));   end   ha = Bus.int(round(SAE1.con(i,1)));   hb = Bus.int(round(SAE1.con(i,2)));   xt1 = SAE1.con(i,6);   xt2 = SAE1.con(i,7);   xt3 = SAE1.con(i,8);   xa0 = SAE1.con(i,9);   xb0 = SAE1.con(i,10);   if xt1 == 0; xt1 = 0.2; end   if xt2 == 0; xt2 = 0.2; end   if xt3 == 0; xt3 = 0.2; end   if xa0 == 0; xa0 = 0.2; end   if xb0 == 0; xb0 = 0.2; end   a1 = SAE1.con(i,11);   b1 = SAE1.con(i,12);   a2 = SAE1.con(i,13);   b2 = SAE1.con(i,14);   a3 = SAE1.con(i,15);   b3 = SAE1.con(i,16);   h1 = SAE1.con(i,17);   k1 = SAE1.con(i,18);   vrif1 = SAE1.con(i,19);   h2 = SAE1.con(i,20);   k2 = SAE1.con(i,21);   vrif2 = SAE1.con(i,22);   h3 = SAE1.con(i,23);   k3 = SAE1.con(i,24);   vrif3 = SAE1.con(i,25);   mmax1 = SAE1.con(i,26);   mmin1 = SAE1.con(i,27);   mmax2 = SAE1.con(i,28);   mmin2 = SAE1.con(i,29);   mmax3 = SAE1.con(i,30);   mmin3 = SAE1.con(i,31);   m0 = (0.9519 + 0.9474 + 0.9492)/3;   mmax = (mmax1+mmax2+mmax3)/3;   mmin = (mmin1+mmin2+mmin3)/3;   %  coefficienti nel caso mi = m + Bi   % B1 = 0.9519 - m0;   % B2 = 0.9474 - m0;   % B3 = 0.9492 - m0;   %  coefficienti nel caso m1 = m per qualunque mi   B1 = 0;   B2 = 0;   B3 = 0;   va = DAE.V(ha);   vb = DAE.V(hb);   delta = DAE.a(ha);   theta = DAE.a(hb);   if isempty(SAE1.m)      hm = 1;   else      hm = SAE1.m(i);   end   t1 = m+B1;   t2 = 1/t1;   t4 = 1/xt1;   t5 = t1*t1;   t6 = 1/t5;   t8 = t4+a1*t6;   t9 = 1/t8;   t10 = t4*t9;   t12 = m+B2;   t13 = 1/t12;   t15 = 1/xt2;   t16 = t12*t12;   t17 = 1/t16;   t19 = t15+a2*t17;   t20 = 1/t19;   t21 = t15*t20;   t23 = m+B3;   t24 = 1/t23;   t26 = 1/xt3;   t27 = t23*t23;   t28 = 1/t27;   t30 = t26+a3*t28;   t31 = 1/t30;   t32 = t26*t31;   t34 = b1*t2*t10+b2*t13*t21+b3*t24*t32;   t35 = 1/xa0;   t36 = 1/xb0;   t37 = xt1*xt1;   t38 = 1/t37;   t40 = xt2*xt2;   t41 = 1/t40;   t43 = xt3*xt3;   t44 = 1/t43;   t46 = t4+t26+t15+t35+t36-t38*t9-t41*t20-t44*t31;   t47 = 1/t46;   t48 = t34*t47;   t49 = va*t35;   t51 = va*vb;   t52 = -delta+theta;   t53 = sin(t52);   t54 = t51*t53;   t55 = t47*t35;   t56 = t55*t36;   t57 = t54*t56;   t59 = vb*t36;   t62 = xa0*xa0;   t63 = 1/t62;   t65 = t35-t63*t47;   t66 = va*va;   t68 = cos(t52);   t70 = t56*t51*t68;   t72 = xb0*xb0;   t73 = 1/t72;   t75 = t36-t73*t47;   t76 = vb*vb;   t81 = vb*t53*t56;   t84 = va*t53*t56;   t87 = t55*t59*t68;   t91 = t55*t36*va*t68;   t97 = t49+t59;   t98 = t97*t47;   t100 = t2*t4*t9;   t106 = t13*t15*t20;   t112 = t24*t26*t31;   t117 = t46*t46;   t118 = 1/t117;   t119 = t97*t118;   t121 = t8*t8;   t122 = 1/t121;   t128 = t19*t19;   t129 = 1/t128;   t135 = t30*t30;   t136 = 1/t135;   t142 = -2.0*t38*t122*a1/t5/t1-2.0*t41*t129*a2/t16/t12-2.0*t44*t136*a3/t27/t23;   t148 = t5*t5;   t149 = 1/t148;   t152 = t4*t122*a1;   t162 = t16*t16;   t163 = 1/t162;   t166 = t15*t129*a2;   t176 = t27*t27;   t177 = 1/t176;   t180 = t26*t136*a3;   t198 = (-b1*t6*t10+2.0*b1*t149*t152-b2*t17*t21+2.0*b2*t163*t166-b3*t28* ...      t32+2.0*b3*t177*t180)*t47;   t200 = t34*t118;   t203 = t118*t35;   t206 = t54*t203*t36*t142;   t214 = t203*t36*t51*t68*t142;   switch flag   case 0      if length(SAE1.con(1,:)) > 28          xeq = zeros(5,1);          A =  [1 0 1 0 0; 1 0 0 1 0; 1 0 0 0 1; 0 1 1 0 0; 0 1 0 1 0; ...                  0 1 0 0 1; 0 0 1 1 0; 0 0 0 1 1; 0 0 1 0 1; 1 1 0 0 0];          x1 = SAE1.con(i,35) + SAE1.con(i,32);          x2 = SAE1.con(i,35) + SAE1.con(i,33) + SAE1.con(i,37);          x3 = SAE1.con(i,35) + SAE1.con(i,34) + SAE1.con(i,37) + SAE1.con(i,38);          x4 = SAE1.con(i,36) + SAE1.con(i,32) + SAE1.con(i,38) + SAE1.con(i,37);          x5 = SAE1.con(i,36) + SAE1.con(i,33) + SAE1.con(i,38);          x6 = SAE1.con(i,36) + SAE1.con(i,34);          x7 = SAE1.con(i,32) + SAE1.con(i,37) + SAE1.con(i,33);          x8 = SAE1.con(i,34) + SAE1.con(i,38) + SAE1.con(i,33);          x9 = SAE1.con(i,32) + SAE1.con(i,37) + SAE1.con(i,38) + SAE1.con(i,34);          x10= SAE1.con(i,35) + SAE1.con(i,37) + SAE1.con(i,38) + SAE1.con(i,36);          xreal = [x1; x2; x3; x4; x5; x6; x7; x8; x9; x10];          xeq = A\xreal;          SAE1.con(i,9) = xeq(1);          SAE1.con(i,10) = xeq(2);          SAE1.con(i,6) = xeq(3);          SAE1.con(i,7) = xeq(4);          SAE1.con(i,8) = xeq(5);      end   case 1      Pa = t48*t49-t57;      Pb = t48*t59+t57;      Qa = t65*t66-t70;      Qb = t75*t76-t70;      DAE.gp(ha) = Pa + DAE.gp(ha);      DAE.gq(ha) = Qa + DAE.gq(ha);      DAE.gp(hb) = Pb + DAE.gp(hb);      DAE.gq(hb) = Qb + DAE.gq(hb);      % xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx   case 2      DAE.J12(ha,ha) = DAE.J12(ha,ha) + t48*t35-t81;      DAE.J12(ha,hb) = DAE.J12(ha,hb) - t84;      DAE.J11(ha,ha) = DAE.J11(ha,ha) + t70;      DAE.J11(ha,hb) = DAE.J11(ha,hb) - t70;      DAE.J22(ha,ha) = DAE.J22(ha,ha) + 2.0*t65*va-t87;      DAE.J22(ha,hb) = DAE.J22(ha,hb) - t91;      DAE.J21(ha,ha) = DAE.J21(ha,ha) - t57;      DAE.J21(ha,hb) = DAE.J21(ha,hb) + t57;      DAE.J12(hb,ha) = DAE.J12(hb,ha) + t81;      DAE.J12(hb,hb) = DAE.J12(hb,hb) + t48*t36+t84;      DAE.J11(hb,ha) = DAE.J11(hb,ha) - t70;      DAE.J11(hb,hb) = DAE.J11(hb,hb) + t70;      DAE.J22(hb,ha) = DAE.J22(hb,ha) - t87;      DAE.J22(hb,hb) = DAE.J22(hb,hb) + 2.0*t75*vb-t91;      DAE.J21(hb,ha) = DAE.J21(hb,ha) -t57;      DAE.J21(hb,hb) = DAE.J21(hb,hb) +  t57;      % xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx   case 3      %   DAE.f :      DAE.f(hm) = -h1*t1/3+k1*(t98*t100-vrif1)/3-h2*t12/3+k2*(t98*t106-vrif2)/3-h3* ...         t23/3+k3*(t98*t112-vrif3)/3;      if (m >= mmax)         m = mmax;         DAE.x(hm) = m;         if(DAE.f(hm) > 0); DAE.f(hm) = 0;end      end      if (m <= mmin)         m = mmin;         DAE.x(hm) = m;         if(DAE.f(hm) < 0); DAE.f(hm) = 0;end      end      % xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx   case 4      DAE.Fx(hm,hm) = -h1/3+k1*(-t119*t2*t10*t142-t98*t6*t4*t9+2.0*t98*t149*t152)/3-h2 ...         /3+k2*(-t119*t13*t21*t142-t98*t17*t15*t20+2.0*t98*t163*t166)/3-h3/3+k3*(-t119* ...         t24*t32*t142-t98*t28*t26*t31+2.0*t98*t177*t180)/3;      DAE.Gx(ha,hm) = t198*t49-t200*t49*t142+t206;      DAE.Gx(ha+Bus.n,hm) = t63*t118*t142*t66+t214;      DAE.Gx(hb,hm) = t198*t59-t200*t59*t142-t206;      DAE.Gx(hb+Bus.n,hm) = t73*t118*t142*t76+t214;      DAE.Fy(hm,Bus.n+ha) = k1*t35*t47*t100/3+k2*t35*t47*t106/3+k3*t35*t47*t112/3;      DAE.Fy(hm,Bus.n+hb) = k1*t36*t47*t100/3+k2*t36*t47*t106/3+k3*t36*t47*t112/3;      if ((m >= mmax | m <= mmin) & DAE.f(hm) == 0)         DAE.Fx(hm,hm) = -1;         DAE.Fy(hm,:) = zeros(1,2*Bus.n);      end      % xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx   case 5      if ((m >= mmax | m <= mmin) & DAE.f(hm) == 0)         DAE.tn(hm) = 0;         DAE.Ac(hm,:) = zeros(1,DAE.n+2*Bus.n);         DAE.Ac(hm,hm) = DAE.Fx(hm,hm);      end   endend