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 Milano

global DAE Bus SAE3

for 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

    end
end