
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%        FORMATION OF A MATRIX with PSS

TAGmod=sparse([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(Ag21)  zeros(ngen)  diag(Ag23)    diag(Ag24)   diag(Ag25)   diag(Ag26)    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)   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)   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)   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)   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)   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)   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)   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)   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)   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)   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)   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)   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)   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)   zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)  zeros(ngen)   zeros(ngen)  zeros(ngen)]);    
              
 TAG=TAGmod(ii,ii);
 clear TAGmod;
 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 TBGmod=sparse([zeros(ngen)      zeros(ngen)
              diag(Bg21rect)   diag(Bg22rect)
              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)
              zeros(ngen)      zeros(ngen)
              zeros(ngen)      zeros(ngen)
              zeros(ngen)      zeros(ngen)
              zeros(ngen)      zeros(ngen)
              zeros(ngen)      zeros(ngen)
              zeros(ngen)      zeros(ngen) ]);
           
           
  TBG=TBGmod(ii,kk);
  clear TBGmod; 
  %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 if(npss~=0)
   TAGp=([TAG+sparse(zeros(size(EG*PS1'*DPS*PS1*FG))) sparse(zeros(size(EG*PS1'*CPS))); BPS*PS1*FG APS]);
   TBGp=([TBG; sparse(zeros(size(APS,1),2*ngen))]);
   TCGp=CG;
   TEGp=([EG; sparse(zeros(size(APS,1),ngen))]); %padded zeros to EG suitable.
   TFGp=([FG  sparse(zeros(ngen,size(APS,1)))]); %padded zeros to FG suitable.
   TEATd=-TAGp-TBGp*PG'*YDQdi*PG*TCGp;
else
   TEGp=EG;
   TFGp=FG; 
  TEATd=-TAG-TBG*PG'*YDQdi*PG*CG;
 end
 %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 H2_diag = 2*diag(H);
 TEATD=H2_diag*TEATd([2:16:16*ngen],:);

APSS=[];
BPSS=[];
CPSS=[];
DPSS=[];

for k=1:size(delPw_pss,1)
[N1,D1]=series([Tw1_delPw(k) 0],[Tw1_delPw(k) 1],[Tw2_delPw(k) 0],[Tw2_delPw(k) 1]);
[N1,D1]=series(N1,D1,[0 1],[T6_delPw(k) 1]);
[APW,BPW,CPW,DPW]=tf2ss(N1,D1);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[N2,D2]=series([Tw3_delPw(k) 0],[Tw3_delPw(k) 1],[Tw4_delPw(k) 0],[Tw4_delPw(k) 1]);
[N2,D2]=series(N2,D2,[0 T7_delPw(k)/2/H_delPw(k)],[T7_delPw(k) 1]);
[APT,BPT,CPT,DPT]=tf2ss(N2,D2);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[N3,D3]=series([T8_delPw(k) 1],[(T9_delPw(k)*T9_delPw(k)) (2*T9_delPw(k)) 1],[T8_delPw(k) 1],[(T9_delPw(k)*T9_delPw(k)) (2*T9_delPw(k)) 1]);
[N3,D3]=series(N3,D3,[T8_delPw(k) 1],[(T9_delPw(k)*T9_delPw(k)) (2*T9_delPw(k)) 1]);
[N3,D3]=series(N3,D3,[T8_delPw(k) 1],[(T9_delPw(k)*T9_delPw(k)) (2*T9_delPw(k)) 1]);
[APF,BPF,CPF,DPF]=tf2ss(N3,D3);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[N4,D4]=series([T1_delPw(k) 1],[T2_delPw(k) 1],[Ks1_delPw(k)],[1]);
[N4,D4]=series(N4,D4,[T3_delPw(k) 1],[T4_delPw(k) 1]);
[APC,BPC,CPC,DPC]=tf2ss(N4,D4);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
APSSd=[APW    zeros(size(APW,1),size(APT,1)) zeros(size(APW,1),size(APF,1))  zeros(size(APW,1),size(APC,1))
    zeros(size(APT,1),size(APW,1))    APT    zeros(size(APT,1),size(APF,1))  zeros(size(APT,1),size(APC,1))
    (BPF*CPW)      (BPF*Ks3_delPw(k)*CPT)                APF                    zeros(size(APF,1),size(APC,1))
    (BPC*DPF*CPW)  (BPC*(DPF*Ks3_delPw(k)*CPT-CPT))    (BPC*CPF)                APC];
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
BPSSd=[     BPW         zeros(size(APW,1),1)
       zeros(size(APT,1),1)     BPT
       (BPF*DPW)       (BPF*Ks3_delPw(k)*DPT)
       (BPC*DPF*DPW)   (BPC*(DPF*Ks3_delPw(k)*DPT-DPT))];
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
CPSSd=[DPC*DPF*CPW   (DPC*(DPF*Ks3_delPw(k)*CPT-CPT))  DPC*CPF    CPC];
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
DPSSd=[DPC*DPF*DPW   (DPC*(DPF*Ks3_delPw(k)*DPT-DPT))];
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~       
[APSS,BPSS,CPSS,DPSS]=append(APSS,BPSS,CPSS,DPSS,APSSd,BPSSd,CPSSd,DPSSd);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
APSS=sparse(APSS);
BPSS=sparse(BPSS);
CPSS=sparse(CPSS);
DPSS=sparse(DPSS);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for k=1:size(delPw_pss,1)
    for p=1:ngen
        if (delPw_pss(k,1)==gen(p,1))
            temp(k)=p;
        end
    end
end

PS2=[];
PS2m=[];
for q=1:size(delPw_pss,1)
PS21=zeros(1,ngen);
PS21(1,temp(q))=1;
PS2m1=[PS21 zeros(1,ngen);zeros(1,ngen) PS21];
PS2=[PS2;PS21];
PS2m=[PS2m;PS2m1];
end

FGd=[TFGp;TEATD];
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
AG=([AG+TEGp*PS2'*DPSS*PS2m*FGd TEGp*PS2'*CPSS; BPSS*PS2m*FGd APSS]);
BG=([BG; sparse(zeros(size(APSS,1),2*ngen))]);
CG=([CG  sparse(zeros(size(APSS,1),2*ngen)')]);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
clear C1 D1;
clear A1 B1;
clear temp;
