function fm_sae3(flag)% FM_SAE3 define a subtransmission area equivalent%         with three loads and three LTCs. This model%         uses three state variables for the tap ratios%         of the LTCs.%% FM_SAE3(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_SAE1 and FM_SAE2%%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 DAE Bus SAE3for i = 1:SAE3.n    if isempty(SAE3.m1)        m1 = 1;        m2 = 1;        m3 = 1;    else        m1 = DAE.x(SAE3.m1(i));        m2 = DAE.x(SAE3.m2(i));        m3 = DAE.x(SAE3.m3(i));    end    ha = Bus.int(round(SAE3.con(i,1)));    hb = Bus.int(round(SAE3.con(i,2)));    xt1 = SAE3.con(i,6);    xt2 = SAE3.con(i,7);    xt3 = SAE3.con(i,8);    xa0 = SAE3.con(i,9);    xb0 = SAE3.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 = SAE3.con(i,11);    b1 = SAE3.con(i,12);    a2 = SAE3.con(i,13);    b2 = SAE3.con(i,14);    a3 = SAE3.con(i,15);    b3 = SAE3.con(i,16);    h1 = SAE3.con(i,17);    k1 = SAE3.con(i,18);    vrif1 = SAE3.con(i,19);    h2 = SAE3.con(i,20);    k2 = SAE3.con(i,21);    vrif2 = SAE3.con(i,22);    h3 = SAE3.con(i,23);    k3 = SAE3.con(i,24);    vrif3 = SAE3.con(i,25);    mmax1 = SAE3.con(i,26);    mmin1 = SAE3.con(i,27);    mmax2 = SAE3.con(i,28);    mmin2 = SAE3.con(i,29);    mmax3 = SAE3.con(i,30);    mmin3 = SAE3.con(i,31);    va = DAE.V(ha);    vb = DAE.V(hb);    delta = DAE.a(ha);    theta = DAE.a(hb);    if isempty(SAE3.m1)        hm1 = 1;        hm2 = 1;        hm3 = 1;    else        hm1 = SAE3.m1(i);        hm2 = SAE3.m2(i);        hm3 = SAE3.m3(i);    end    t1 = 1/m1;    t3 = 1/xt1;    t4 = m1*m1;    t5 = 1/t4;    t7 = t3+a1*t5;    t8 = 1/t7;    t9 = t3*t8;    t11 = 1/m2;    t13 = 1/xt2;    t14 = m2*m2;    t15 = 1/t14;    t17 = t13+a2*t15;    t18 = 1/t17;    t19 = t13*t18;    t21 = 1/m3;    t23 = 1/xt3;    t24 = m3*m3;    t25 = 1/t24;    t27 = t23+a3*t25;    t28 = 1/t27;    t29 = t23*t28;    t31 = b1*t1*t9+b2*t11*t19+b3*t21*t29;    t32 = 1/xa0;    t33 = 1/xb0;    t34 = xt1*xt1;    t35 = 1/t34;    t37 = xt2*xt2;    t38 = 1/t37;    t40 = xt3*xt3;    t41 = 1/t40;    t43 = t3+t23+t13+t32+t33-t35*t8-t38*t18-t41*t28;    t44 = 1/t43;    t45 = t31*t44;    t46 = va*t32;    t48 = va*vb;    t49 = delta-theta;    t50 = sin(t49);    t52 = t44*t32;    t53 = t52*t33;    t54 = t48*t50*t53;    t56 = vb*t33;    t59 = xa0*xa0;    t60 = 1/t59;    t62 = t32-t60*t44;    t63 = va*va;    t65 = cos(t49);    t67 = t53*t48*t65;    t69 = xb0*xb0;    t70 = 1/t69;    t72 = t33-t70*t44;    t73 = vb*vb;    t78 = vb*t50*t53;    t81 = va*t50*t53;    t84 = t52*t56*t65;    t86 = t33*va;    t88 = t52*t86*t65;    t94 = t46+t56;    t95 = t94*t44;    t97 = t1*t3*t8;    t104 = t11*t13*t18;    t111 = t21*t23*t28;    t116 = t43*t43;    t117 = 1/t116;    t118 = t94*t117;    t119 = t4*t4;    t120 = 1/t119;    t124 = t7*t7;    t134 = 1/t124;    t136 = t3*t134*a1;    t144 = k1*t94*t117*t1*t3;    t146 = t17*t17;    t147 = 1/t146;    t148 = t147*a2;    t150 = 1/t14/m2;    t151 = t148*t150;    t155 = t27*t27;    t156 = 1/t155;    t157 = t156*a3;    t159 = 1/t24/m3;    t160 = t157*t159;    t166 = k2*t94*t117*t11*t13;    t168 = t134*a1;    t170 = 1/t4/m1;    t171 = t168*t170;    t174 = t14*t14;    t175 = 1/t174;    t189 = t13*t147*a2;    t200 = k3*t94*t117*t21*t23;    t207 = t24*t24;    t208 = 1/t207;    t222 = t23*t156*a3;    t232 = (-b1*t5*t9+2.0*b1*t120*t136)*t44;    t234 = t31*t117;    t235 = t234*t46;    t238 = t35*t134*a1*t170;    t242 = t48*t50*t117*t32;    t245 = t242*t33*t35*t171;    t252 = (-b2*t15*t19+2.0*b2*t175*t189)*t44;    t256 = t38*t147*a2*t150;    t260 = t242*t33*t38*t151;    t267 = (-b3*t25*t29+2.0*b3*t208*t222)*t44;    t271 = t41*t156*a3*t159;    t275 = t242*t33*t41*t160;    t277 = t60*t117;    t284 = t32*t117*t86*vb;    t287 = t284*t65*t35*t171;    t295 = t284*t65*t38*t151;    t303 = t284*t65*t41*t160;    t306 = t234*t56;    t315 = t70*t117;    switch flag    case 0        if length(SAE3.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 = SAE3.con(i,35) + SAE3.con(i,32);            x2 = SAE3.con(i,35) + SAE3.con(i,33) + SAE3.con(i,37);            x3 = SAE3.con(i,35) + SAE3.con(i,34) + SAE3.con(i,37) + SAE3.con(i,38);            x4 = SAE3.con(i,36) + SAE3.con(i,32) + SAE3.con(i,38) + SAE3.con(i,37);            x5 = SAE3.con(i,36) + SAE3.con(i,33) + SAE3.con(i,38);            x6 = SAE3.con(i,36) + SAE3.con(i,34);            x7 = SAE3.con(i,32) + SAE3.con(i,37) + SAE3.con(i,33);            x8 = SAE3.con(i,34) + SAE3.con(i,38) + SAE3.con(i,33);            x9 = SAE3.con(i,32) + SAE3.con(i,37) + SAE3.con(i,38) + SAE3.con(i,34);            x10= SAE3.con(i,35) + SAE3.con(i,37) + SAE3.con(i,38) + SAE3.con(i,36);            xreal = [x1; x2; x3; x4; x5; x6; x7; x8; x9; x10];            xeq = A\xreal;            SAE3.con(i,9) = xeq(1);            SAE3.con(i,10) = xeq(2);            SAE3.con(i,6) = xeq(3);            SAE3.con(i,7) = xeq(4);            SAE3.con(i,8) = xeq(5);        end    case 1        Pa = t45*t46+t54;        Pb = t45*t56-t54;        Qa = t62*t63-t67;        Qb = t72*t73-t67;        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);    case 2        DAE.J12(ha,ha) = DAE.J12(ha,ha) + t45*t32+t78;        DAE.J12(ha,hb) = DAE.J12(ha,hb) + t81;        DAE.J11(ha,ha) = DAE.J11(ha,ha) + t67;        DAE.J11(ha,hb) = DAE.J11(ha,hb) - t67;        DAE.J22(ha,ha) = DAE.J22(ha,ha) + 2.0*t62*va-t84;        DAE.J22(ha,hb) = DAE.J22(ha,hb) - t88;        DAE.J21(ha,ha) = DAE.J21(ha,ha) + t54;        DAE.J21(ha,hb) = DAE.J21(ha,hb) - t54;        DAE.J12(hb,ha) = DAE.J12(hb,ha) - t78;        DAE.J12(hb,hb) = DAE.J12(hb,hb) + t45*t33-t81;        DAE.J11(hb,ha) = DAE.J11(hb,ha) - t67;        DAE.J11(hb,hb) = DAE.J11(hb,hb) + t67;        DAE.J22(hb,ha) = DAE.J22(hb,ha) - t84;        DAE.J22(hb,hb) = DAE.J22(hb,hb) + 2.0*t72*vb-t88;        DAE.J21(hb,ha) = DAE.J21(hb,ha) + t54;        DAE.J21(hb,hb) = DAE.J21(hb,hb) - t54;    case 3        DAE.f(hm1) = -h1*m1+k1*(t95*t97-vrif1);        DAE.f(hm2) = -h2*m2+k2*(t95*t104-vrif2);        DAE.f(hm3) = -h3*m3+k3*(t95*t111-vrif3);        if (m1 >= mmax1)            m1 = mmax1;            DAE.x(hm1) = m1;            if(DAE.f(hm1) > 0); DAE.f(hm1) = 0;end        end        if (m1 <= mmin1)            m1 = mmin1;            DAE.x(hm1) = m1;            if(DAE.f(hm1) < 0); DAE.f(hm1) = 0;end        end        if (m2 >= mmax2)            m2 = mmax2;            DAE.x(hm2) = m2;            if(DAE.f(hm2) > 0); DAE.f(hm2) = 0;end        end        if (m2 <= mmin2)            m2 = mmin2;            DAE.x(hm2) = m2;            if(DAE.f(hm2) < 0); DAE.f(hm2) = 0;end        end        if (m3 >= mmax3)            m3 = mmax3;            DAE.x(hm3) = m3;            if(DAE.f(hm3) > 0); DAE.f(hm3) = 0;end        end        if (m3 <= mmin3)            m3 = mmin3;            DAE.x(hm3) = m3;            if(DAE.f(hm3) < 0); DAE.f(hm3) = 0;end        end    case 4        DAE.Fx(hm1,hm1) = -h1+k1*(2.0*t118*t120/t34/xt1/t124/t7*a1-t95*t5*t3*t8+2.0*t95*t120*t136);        DAE.Fx(hm1,hm2) = 2.0*t144*t8*t38*t151;        DAE.Fx(hm1,hm3) = 2.0*t144*t8*t41*t160;        DAE.Fx(hm2,hm1) = 2.0*t166*t18*t35*t171;        DAE.Fx(hm2,hm2) = -h2+k2*(2.0*t118*t175/t37/xt2/t146/t17*a2-t95*t15*t13*t18+2.0*t95*t175*t189);        DAE.Fx(hm2,hm3) = 2.0*t166*t18*t41*t160;        DAE.Fx(hm3,hm1) = 2.0*t200*t28*t35*t171;        DAE.Fx(hm3,hm2) = 2.0*t200*t28*t38*t151;        DAE.Fx(hm3,hm3) = -h3+k3*(2.0*t118*t208/t40/xt3/t155/t27*a3-t95*t25*t23*t28+2.0*t95*t208*t222);        DAE.Gx(ha,hm1) = t232*t46+2.0*t235*t238+2.0*t245;        DAE.Gx(ha,hm2) = t252*t46+2.0*t235*t256+2.0*t260;        DAE.Gx(ha,hm3) = t267*t46+2.0*t235*t271+2.0*t275;        DAE.Gx(ha+Bus.n,hm1) = -2.0*t277*t35*t168*t170*t63-2.0*t287;        DAE.Gx(ha+Bus.n,hm2) = -2.0*t277*t38*t148*t150*t63-2.0*t295;        DAE.Gx(ha+Bus.n,hm3) = -2.0*t277*t41*t157*t159*t63-2.0*t303;        DAE.Gx(hb,hm1) = t232*t56+2.0*t306*t238-2.0*t245;        DAE.Gx(hb,hm2) = t252*t56+2.0*t306*t256-2.0*t260;        DAE.Gx(hb,hm3) = t267*t56+2.0*t306*t271-2.0*t275;        DAE.Gx(hb+Bus.n,hm1) = -2.0*t315*t35*t168*t170*t73-2.0*t287;        DAE.Gx(hb+Bus.n,hm2) = -2.0*t315*t38*t148*t150*t73-2.0*t295;        DAE.Gx(hb+Bus.n,hm3) = -2.0*t315*t41*t157*t159*t73-2.0*t303;        DAE.Fy(hm1,Bus.n+ha) = k1*t32*t44*t97;        DAE.Fy(hm1,Bus.n+hb) = k1*t33*t44*t97;        DAE.Fy(hm2,Bus.n+ha) = k2*t32*t44*t104;        DAE.Fy(hm2,Bus.n+hb) = k2*t33*t44*t104;        DAE.Fy(hm3,Bus.n+ha) = k3*t32*t44*t111;        DAE.Fy(hm3,Bus.n+hb) = k3*t33*t44*t111;        if ((m1 >= mmax1 | m1 <= mmin1) & DAE.f(hm1) == 0)            DAE.Fx(hm1,:) = zeros(1,DAE.n);            DAE.Fx(hm1,hm1) = -1;            DAE.Fy(hm1,:) = zeros(1,2*Bus.n);        end        if ((m2 >= mmax2 | m2 <= mmin2) & DAE.f(hm2) == 0)            DAE.Fx(hm2,:) = zeros(1,DAE.n);            DAE.Fx(hm2,hm2) = -1;            DAE.Fy(hm2,:) = zeros(1,2*Bus.n);        end        if ((m3 >= mmax3 | m3 <= mmin3) & DAE.f(hm3) == 0)            DAE.Fx(hm3,:) = zeros(1,DAE.n);            DAE.Fx(hm3,hm3) = -1;            DAE.Fy(hm3,:) = zeros(1,2*Bus.n);        end    case 5        if ((m1 >= mmax1 | m1 <= mmin1) & DAE.f(hm1) == 0)            DAE.tn(hm1) = 0;            DAE.Ac(hm1,:) = zeros(1,DAE.n+2*Bus.n);            DAE.Ac(hm1,hm1) = DAE.Fx(hm1,hm1);        end        if ((m2 >= mmax2 | m2 <= mmin2) & DAE.f(hm2) == 0)            DAE.tn(hm2) = 0;            DAE.Ac(hm2,:) = zeros(1,DAE.n+2*Bus.n);            DAE.Ac(hm2,hm2) = DAE.Fx(hm2,hm2);        end        if ((m3 >= mmax3 | m3 <= mmin3) & DAE.f(hm3) == 0)            DAE.tn(hm3) = 0;            DAE.Ac(hm3,:) = zeros(1,DAE.n+2*Bus.n);            DAE.Ac(hm3,hm3) = DAE.Fx(hm3,hm3);        end    endend