
% Formation of the generator matrices


%****************************************************************************************************
%****************************************************************************************************
%             FORMATION OF THE `A' MATRIX BEGINS
%****************************************************************************************************
%****************************************************************************************************

dmth0=delta0-thetag0;
C1=(xdd-xddd)./xdd;
C2=(xd-xdd).*xddd./xd./xdd;
C3=(xqd-xqdd)./xqd;
C4=(xq-xqd).*xqdd./xq./xqd;
C5=(xddd-xqdd)./xddd./xqdd;

Ag12=wB*ones(1,ngen);

Ag21=(C2.*Vg0.*sif0.*cos(dmth0)./xddd + C1.*Vg0.*sih0.*cos(dmth0)./xddd +C4.*Vg0.*sig0.*sin(dmth0)./xqdd);
Ag21=(Ag21+C3.*Vg0.*sik0.*sin(dmth0)./xqdd + C5.*Vg0.*Vg0.*cos(2.0*(dmth0)))*(-0.5)./H;
Ag22=-(Dm./H)*(0.5);

Ag23=-C2.*Vg0.*sin(dmth0)*(0.5)./H./xddd;

Ag24=-C1.*Vg0.*sin(dmth0)*(0.5)./H./xddd;

Ag25=C4.*Vg0.*cos(dmth0)*(0.5)./H./xqdd;

Ag26=C3.*Vg0.*cos(dmth0)*(0.5)./H./xqdd;

Ag31=-Vg0.*sin(dmth0)./Tdd;

Ag33=-ones(1,ngen)./Tdd;

Ag37=(xdd./(xd-xdd))./Tdd;

Ag41=-Vg0.*sin(dmth0)./Tddd;

Ag44=-ones(1,ngen)./Tddd;

Ag51=Vg0.*cos(dmth0)./Tqd;

Ag55=-ones(1,ngen)./Tqd;

Ag61=Vg0.*cos(dmth0)./Tqdd;

Ag66=-ones(1,ngen)./Tqdd;

Ag77=static_dEfd+DC1A_dAg77+AC4A_dAg77;

Ag78=DC1A_dAg78;

Ag79=AC4A_dAg79;

Ag87=DC1A_dAg87;

Ag88=DC1A_dAg88;

Ag89=DC1A_dAg89;

Ag810=DC1A_dAg810;

Ag97=DC1A_dAg97;

Ag99=DC1A_dAg99+AC4A_dAg99;

Ag910=DC1A_dAg910;

Ag107=DC1A_dAg107;

Ag1010=DC1A_dAg1010;

Ag211=RHST_dAg211;

Ag212=RHST_dAg212;

Ag213=RHST_dAg213;

Ag215=HYDRO_dAg215;

Ag216=HYDRO_dAg216;

Ag1111=RHST_dAg1111;

Ag1115=RHST_dAg1115;

Ag1211=RHST_dAg1211;

Ag1212=RHST_dAg1212;

Ag1312=RHST_dAg1312;

Ag1313=RHST_dAg1313;

Ag142=RHST_dAg142+HYDRO_dAg142;

Ag1414=RHST_dAg1414+HYDRO_dAg1414;

Ag152=RHST_dAg152+HYDRO_dAg152;

Ag1514=RHST_dAg1514+HYDRO_dAg1514;

Ag1515=RHST_dAg1515+HYDRO_dAg1515;

Ag1615=HYDRO_dAg1615;

Ag1616=HYDRO_dAg1616;




AGmod=sparse([zeros(ngen) diag(Ag12)  zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen)  zeros(ngen)
              diag(Ag21)  diag(Ag22)  diag(Ag23)    diag(Ag24)   diag(Ag25)   diag(Ag26)    zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  diag(Ag211)  diag(Ag212)  diag(Ag213)  zeros(ngen)   diag(Ag215)  diag(Ag216)
              diag(Ag31)  zeros(ngen) diag(Ag33)    zeros(ngen)  zeros(ngen)  zeros(ngen)   diag(Ag37)  zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen)  zeros(ngen)
              diag(Ag41)  zeros(ngen) zeros(ngen)   diag(Ag44)   zeros(ngen)  zeros(ngen)   zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen)  zeros(ngen)
              diag(Ag51)  zeros(ngen) zeros(ngen)   zeros(ngen)  diag(Ag55)   zeros(ngen)   zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen)  zeros(ngen)
              diag(Ag61)  zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  diag(Ag66)    zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen)  zeros(ngen)
              zeros(ngen) zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)   diag(Ag77)  diag(Ag78)    diag(Ag79)  zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen)  zeros(ngen)
              zeros(ngen) zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)   diag(Ag87)  diag(Ag88)    diag(Ag89)   diag(Ag810)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen)  zeros(ngen)
              zeros(ngen) zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)   diag(Ag97)  zeros(ngen)   diag(Ag99)   diag(Ag910)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen)  zeros(ngen)
              zeros(ngen) zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)   diag(Ag107) zeros(ngen)   zeros(ngen)  diag(Ag1010) zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen)  zeros(ngen)
              zeros(ngen) zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  diag(Ag1111) zeros(ngen)  zeros(ngen)  zeros(ngen)   diag(Ag1115) zeros(ngen)                                        
              zeros(ngen) zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  diag(Ag1211) diag(Ag1212) zeros(ngen)  zeros(ngen)   zeros(ngen)  zeros(ngen)
              zeros(ngen) zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)   zeros(ngen) diag(Ag1312) diag(Ag1313) zeros(ngen)   zeros(ngen)  zeros(ngen)
              zeros(ngen) diag(Ag142) zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)   zeros(ngen) zeros(ngen)  zeros(ngen)  diag(Ag1414)  zeros(ngen)  zeros(ngen)
              zeros(ngen) diag(Ag152) zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)   zeros(ngen) zeros(ngen)  zeros(ngen)  diag(Ag1514)  diag(Ag1515) zeros(ngen)
              zeros(ngen) zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen) zeros(ngen)   zeros(ngen)  zeros(ngen)   zeros(ngen) zeros(ngen)  zeros(ngen)  zeros(ngen)   diag(Ag1615) diag(Ag1616)]); 
             
           
           
   i=[1:ngen
  (ngen+1):(2*ngen)
  (2*ngen+1):(3*ngen)
  (3*ngen+1):(4*ngen)       
  (4*ngen+1):(5*ngen)
  (5*ngen+1):(6*ngen)
  (6*ngen+1):(7*ngen)
  (7*ngen+1):(8*ngen)
  (8*ngen+1):(9*ngen)
  (9*ngen+1):(10*ngen)
  (10*ngen+1):(11*ngen)
  (11*ngen+1):(12*ngen)
  (12*ngen+1):(13*ngen)
  (13*ngen+1):(14*ngen)
  (14*ngen+1):(15*ngen)
  (15*ngen+1):(16*ngen)];

ii=i(:);

AG=AGmod(ii,ii);

clear AGmod;
%clear Ag21 Ag22 Ag23 Ag24 Ag25 Ag26
clear Ag12  Ag31 Ag33 Ag37 Ag41 Ag44  Ag51 Ag55 Ag61 Ag66 Ag77 Ag78 Ag87 Ag88 Ag89 Ag810 Ag97 Ag99 Ag910 Ag107 Ag1010;
clear Ag211 Ag212 Ag213 Ag214 Ag215 Ag216 Ag1111 Ag1115 Ag1211 Ag1212 Ag1312 Ag1313 Ag142 Ag1414 Ag152 Ag1514 Ag1515 Ag1615 Ag1616; 
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Bg21=(C2.*sif0.*cos(dmth0)./xddd + C1.*sih0.*cos(dmth0)./xddd +C4.*sig0.*sin(dmth0)./xqdd);

Bg21=(Bg21+C3.*sik0.*sin(dmth0)./xqdd + C5.*Vg0.*cos(2.0*(dmth0)))*(0.5)./H;

Bg22=(C2.*sif0.*sin(dmth0)./xddd + C1.*sih0.*sin(dmth0)./xddd -C4.*sig0.*cos(dmth0)./xqdd);

Bg22=(Bg22-C3.*sik0.*cos(dmth0)./xqdd + C5.*Vg0.*sin(2.0*(dmth0)))*(-0.5)./H;

Bg31=sin(dmth0)./Tdd;

Bg32=cos(dmth0)./Tdd;

Bg41=sin(dmth0)./Tddd;

Bg42=cos(dmth0)./Tddd;

Bg51=-cos(dmth0)./Tqd;

Bg52=sin(dmth0)./Tqd;

Bg61=-cos(dmth0)./Tqdd;

Bg62=sin(dmth0)./Tqdd;

Bg72=-static_bdEfd+AC4A_dBg72;

Bg82=DC1A_dBg82;

Bg92=DC1A_dBg92+AC4A_dBg92;

Bg21rect=(Bg21.*(-VDg0)+Bg22.*VQg0)./Vg0;
Bg22rect=(Bg21.*VQg0+Bg22.*VDg0)./Vg0;

Bg31rect=(Bg31.*(-VDg0)+Bg32.*VQg0)./Vg0;
Bg32rect=(Bg31.*VQg0+Bg32.*VDg0)./Vg0;

Bg41rect=(Bg41.*(-VDg0)+Bg42.*VQg0)./Vg0;
Bg42rect=(Bg41.*VQg0+Bg42.*VDg0)./Vg0;

Bg51rect=(Bg51.*(-VDg0)+Bg52.*VQg0)./Vg0;
Bg52rect=(Bg51.*VQg0+Bg52.*VDg0)./Vg0;

Bg61rect=(Bg61.*(-VDg0)+Bg62.*VQg0)./Vg0;
Bg62rect=(Bg61.*VQg0+Bg62.*VDg0)./Vg0;

Bg71rect=(Bg72.*VQg0)./Vg0;
Bg72rect=(Bg72.*VDg0)./Vg0;


Bg81rect=(Bg82.*VQg0)./Vg0;
Bg82rect=(Bg82.*VDg0)./Vg0;


Bg91rect=(Bg92.*VQg0)./Vg0;
Bg92rect=(Bg92.*VDg0)./Vg0;



BGmod=sparse([zeros(ngen)      zeros(ngen)
              diag(Bg21rect)   diag(Bg22rect)
              diag(Bg31rect)   diag(Bg32rect)
              diag(Bg41rect)   diag(Bg42rect)
              diag(Bg51rect)   diag(Bg52rect)
              diag(Bg61rect)   diag(Bg62rect)
              diag(Bg71rect)   diag(Bg72rect)
              diag(Bg81rect)   diag(Bg82rect)
              diag(Bg91rect)   diag(Bg92rect)
              zeros(ngen)      zeros(ngen)
              zeros(ngen)      zeros(ngen)
              zeros(ngen)      zeros(ngen)
              zeros(ngen)      zeros(ngen)
              zeros(ngen)      zeros(ngen)
              zeros(ngen)      zeros(ngen)
              zeros(ngen)      zeros(ngen) ]);
        
       

BG=BGmod(ii,kk);

clear BGmod;
%clear Bg21 Bg22
clear  Bg31 Bg32 Bg41 Bg42 Bg51 Bg52 Bg61 Bg62 Bg72 Bg82 Bg92;
%clear Bg21rect Bg22rect
clear  Bg31rect Bg32rect Bg41rect Bg42rect Bg51rect Bg52rect Bg61rect Bg62rect Bg72rect;
clear Bg81rect Bg82rect Bg91rect Bg92rect;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%Input matrix used to interface PSS -see pssmat.m
EG=sparse(zeros(ngen*16,ngen));
EG([7:16:ngen*16],[1:ngen])=diag(static_bdEfd+AC4A_dEg7Efd);
EG([8:16:ngen*16],[1:ngen])=diag(DC1A_dEg8vr);
EG([9:16:ngen*16],[1:ngen])=diag(DC1A_dEg9xb+AC4A_dEg9xb);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%This matrix is useful to extract 'del s_m' from the 
%x=[deldelta,dels_m,delSif,delSih,delSig,delSik,delEfd,delvr,delxb,delxf,delx1,delx2,delx3,dely1,delPgv,delz1 ]'
FG=sparse(zeros(ngen,ngen*16));
FG([1:ngen],(16*[1:ngen]-14))=eye(ngen);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~       
Cg11=iQg0+Vg0.*(sin(delta0).*cos(dmth0)./xqdd-sin(dmth0).*cos(delta0)./xddd);    
 
Cg13=-C2.*cos(delta0)./xddd;

Cg14=-C1.*cos(delta0)./xddd;

Cg15=-C4.*sin(delta0)./xqdd;

Cg16=-C3.*sin(delta0)./xqdd;

Cg21=Vg0.*(sin(delta0).*sin(dmth0)./xddd+cos(delta0).*cos(dmth0)./xqdd)-iDg0; 
     
Cg23=C2.*sin(delta0)./xddd;

Cg24=C1.*sin(delta0)./xddd;

Cg25=-C4.*cos(delta0)./xqdd;

Cg26=-C3.*cos(delta0)./xqdd;


CGmod=sparse([diag(Cg11) zeros(ngen) diag(Cg13) diag(Cg14)  diag(Cg15)  diag(Cg16)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen) zeros(ngen) zeros(ngen) zeros(ngen) zeros(ngen) zeros(ngen) zeros(ngen)
              diag(Cg21) zeros(ngen) diag(Cg23) diag(Cg24)  diag(Cg25)  diag(Cg26)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen) zeros(ngen) zeros(ngen) zeros(ngen) zeros(ngen) zeros(ngen) zeros(ngen)]);
       
CG=CGmod(kk,ii);

      clear CGmod;
clear Cg11 Cg21 Cg13 Cg14 Cg15 Cg16 Cg23 Cg24  Cg25 Cg26;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


Dg11=cos(delta0).*sin(dmth0)./xddd-sin(delta0).*cos(dmth0)./xqdd;

Dg12=cos(delta0).*cos(dmth0)./xddd+sin(delta0).*sin(dmth0)./xqdd;

Dg21=-sin(delta0).*sin(dmth0)./xddd-cos(delta0).*cos(dmth0)./xqdd;

Dg22=-sin(delta0).*cos(dmth0)./xddd+cos(delta0).*sin(dmth0)./xqdd;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Dg11rect=(Dg11.*(-VDg0)+Dg12.*VQg0)./Vg0;
Dg12rect=(Dg11.*VQg0+Dg12.*VDg0)./Vg0;
Dg21rect=(Dg21.*(-VDg0)+Dg22.*VQg0)./Vg0;
Dg22rect=(Dg21.*VQg0+Dg22.*VDg0)./Vg0;

DGmod=sparse([diag(Dg11rect) diag(Dg12rect); diag(Dg21rect) diag(Dg22rect)]);
 
DG=DGmod(kk,kk);
YG=-DG;
clear DGmod DG;
clear Dg11rect Dg22rect Dg21rect Dg12rect;
clear Dg11 Dg22 Dg21 Dg12;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
pack;

% Now back to the main program
