%pss_selection;
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)  diag(Ag22)  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; 
  %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  TEATd=-TAG-TBG*PG'*YDQdi*PG*CGnp;
  TEAT=TEATd([2:16:16*ngen],:);
  TEAT(:,[1:16:16*ngen])=[]; % eliminating delta for all machines
  TEAT(:,[1:15:15*ngen])=[]; % eliminating Sm for all machines
  rr=input('Enter the generator number for which you want to obtain the angle of GEPS(s):  ');
  for gg=1:1:ngen
     if(gen(gg,1)==rr)
        te=gg;
     end
  end
  TEATR=TEAT(te,:)*2*H(te);%taking the required machine Te data.
  
  %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  TA=Anp; % Anp represents the overall state matrix is obtained in the small_sig programme without PSS. 
  
  TA([1:16:16*ngen],:)=[]; % eliminating delta (in row-wise) for all machines 
  TA([1:15:15*ngen],:)=[]; % eliminating Sm (in row-wise)for all machines
  TA(:,[1:16:16*ngen])=[]; % eliminating delta (in column-wise)for all machines
  TA(:,[1:15:15*ngen])=[]; % eliminating Sm (in column-wise)for all machines
  
  %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  EGT=EG(:,te);      % EG formed in the main programme is used.
  EGT([1:16:16*ngen],:)=[]; % eliminating delta for all machines
  EGT([1:15:15*ngen],:)=[]; % eliminating Sm for all machines
  %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  for i=1:470
    w(i)=0.05*i;
    Xwd(:,i)=(j*w(i)*eye((16-2)*ngen)-TA)\EGT;
  end
  Te_Vs=TEATR*Xwd;
  
  PHASE1=angle(Te_Vs)*180/pi;
  %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 disp('Enter 1  :To design single input PSS - Slip signal ')
 disp('      2  :To design double input PSS ')
 disp('      3  :To design single input PSS -Power signal ')
 psstype = input('Enter your choice : ');
 if(psstype==1)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
mpssd=1;
while(mpssd==1)
 delete(get(0,'children'))
    figure(1)
 plot( w/2/pi, PHASE1);
W_m=input('Enter the center frequency f_m for the PSS (in Hz) :  ' );
Phi_m=input('Enter the amount of phase lead required (in Degrees):  ');
Kslip_pssf=input('Enter the PSS gain Ks:  ' );
alpha=(1-sin(Phi_m*pi/180))/(1+sin(Phi_m*pi/180));
T1_slipf=1/(sqrt(alpha)*2*pi*W_m);
T2_slipf=alpha*T1_slipf;
if(1/alpha<10)
    fprintf('\n The ratio of T1 to T2 = %5.4f is less than 10\n', 1/alpha);
else
    fprintf('\n The ratio of T1 to T2 = %5.4f is greater than 10\n', 1/alpha);
end
 disp('Enter 1 : Only compensator')
 disp('      2 : Washout only')
 disp('      3 : Washout and measuring ckt.')
 disp('      4 : Washout, measuring ckt. and torsional filter')
 config = input('Enter your choice : ');
 
  if (config==1) 
    TFPSS=Kslip_pssf*((1+(j*w).*T1_slipf)./(1+(j*w).*T2_slipf));
  end

  if (config==2)
     Tw_slipf=input('Enter the value of Tw (in s) for the wash-out circuit [1 - 20]s :  ' );
     TFPSS=Kslip_pssf*(((j*w).*Tw_slipf)./(1+(j*w).*Tw_slipf)).*((1+(j*w).*T1_slipf)./(1+(j*w).*T2_slipf));
  end
   
  if (config==3)
   Tw_slipf=input('Enter the value of Tw (in s) for the wash-out circuit [1 - 20]s :  ' );
   TR_slipf=input('Enter the value of TR (in s) for the measurement circuit [0.02-0.05]s :  ' );
   TFPSS=Kslip_pssf*(1./(1+(j*w).*TR_slipf)).*(((j*w).*Tw_slipf)./(1+(j*w).*Tw_slipf)).*((1+(j*w).*T1_slipf)./(1+(j*w).*T2_slipf));
  end
 
  if (config==4)
   Tw_slipf=input('Enter the value of Tw (in s) for the wash-out circuit [1 - 20]s :  ' );
   TR_slipf=input('Enter the value of TR (in s) for the measurement circuit [0.02-0.05]s :  ' );  
   a0 = input('Enter the value of a0 [570] :');
   a1 = input('Enter the value of a1 [35] :');
   TFPSS=Kslip_pssf*(1./(1+(j*w).*TR_slipf)).*(((j*w).*Tw_slipf)./(1+(j*w).*Tw_slipf)).*((1+(j*w).*T1_slipf)./(1+(j*w).*T2_slipf)).*(a0./(((j*w).*(j*w))+a1.*(j*w)+570));
  end 
 
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
PHASE2=angle(TFPSS)*180/pi;
PHASEC=PHASE1+PHASE2;
figure(2)
plot(w/2/pi, PHASEC, 'r', w/2/pi, PHASE1, 'b', w/2/pi, PHASE2, 'k');
grid on
clear PHASEC  PHASE2;
mpssd=input('Enter 0 if you are satisfied with the design, otherwise 1 to re-design the pss:  ');
%-------------------------------------------------------------------

end%while loop end if you are satisfied...psstype==1
if (config==2)
fprintf('\n Update the slip_pss.dat for the machine %d \n',rr);
disp('---------------------------------------------------------')
disp('gen.no   Ks      Tw        T1             T2 ')
disp('---------------------------------------------------------')
fprintf('%d \t %d \t %d \t  %5.5f \t %5.5f \n\n', rr, Kslip_pssf, Tw_slipf, T1_slipf, T2_slipf)
disp('---------------------------------------------------------')
disp('Use typical value for TR (0.02 s), a0 =570, and a1 = 35')
disp('---------------------------------------------------------')
end

if (config==3)
fprintf('\n Update the slip_pss.dat for the machine %d \n',rr);
disp('---------------------------------------------------------')
disp('gen.no   TR              Ks      Tw        T1             T2 ')
disp('---------------------------------------------------------')
fprintf('%d \t %5.5f \t %d \t %d \t  %5.5f \t %5.5f \n\n', rr, TR_slipf, Kslip_pssf, Tw_slipf, T1_slipf, T2_slipf)
disp('---------------------------------------------------------')
disp('Use typical value for a0 =570, and a1 = 35')
disp('---------------------------------------------------------')
end

if (config==4)
fprintf('\n Update the slip_pss.dat for the machine %d \n',rr);
disp('-------------------------------------------------------------------------')
disp('gen.no   TR      Ks      Tw        T1         T2       a0       a1')
disp('-------------------------------------------------------------------------')
fprintf('%d \t %5.3f \t %d \t %d \t  %8.5f  %8.5f  %5.2f  %5.2f \n\n', rr, TR_slipf, Kslip_pssf, Tw_slipf, T1_slipf, T2_slipf,a0,a1)
disp('-------------------------------------------------------------------------')
end

end %if(psstype==1) loop end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if(psstype==2)   
  mpssd=1;
while(mpssd==1)
 delete(get(0,'children'))
 figure(1)
 plot( w/2/pi, PHASE1);
 disp('Now you are designing compensator Gc1(s)..............')
 W_m=input('Enter the center frequency f_m for  Gc1 (in Hz) :  ' );
 Phi_m=input('Enter the amount of phase lead required (in Degrees):  ');
 KdelPw_pssf=input('Enter the PSS gain Ks:  ' );
 alpha=(1-sin(Phi_m*pi/180))/(1+sin(Phi_m*pi/180));
 T1_delPwf=1/(sqrt(alpha)*2*pi*W_m);
 T2_delPwf=alpha*T1_delPwf;
if(1/alpha<10)
    fprintf('\n The ratio of T1 to T2 = %5.4f is less than 10\n', 1/alpha);
else
    fprintf('\n The ratio of T1 to T2 = %5.4f is greater than 10\n', 1/alpha);
 end
 
disp('Enter y : if compensator Gc2(s) is same as the compensator Gc1(s)')
disp('      n : if compensator Gc2(s) is different from the compensator Gc1(s) ')
compensat3= input('Enter your choice (as a character input y or n): '); 
if ( compensat3=='y') 
 T3_delPwf=T1_delPwf;
 T4_delPwf=T2_delPwf;
else
 W_m=input('Enter the center frequency f_m for Gc2 (in Hz) :  ' );
 Phi_m=input('Enter the amount of phase lead required (in Degrees):  ');
 alpha=(1-sin(Phi_m*pi/180))/(1+sin(Phi_m*pi/180));
 T3_delPwf=1/(sqrt(alpha)*2*pi*W_m);
 T4_delPwf=alpha*T3_delPwf;
if(1/alpha<10)
    fprintf('\n The ratio of T3 to T4 = %5.4f is less than 10\n', 1/alpha);
else
    fprintf('\n The ratio of T3 to T4 = %5.4f is greater than 10\n', 1/alpha);
 end
end % if( compensat3=='y') loop end

 disp('Enter 1 : Compensators Gc1(s)*Gc2(s) only')
 disp('      2 : Compensators with Washout only')
 disp('      3 : All Blocks')
 config2 = input('Enter your choice : ');
 if (config2==1) 
    TFPSS=KdelPw_pssf*((1+(j*w).*T1_delPwf)./(1+(j*w).*T2_delPwf)).*((1+(j*w).*T3_delPwf)./(1+(j*w).*T4_delPwf));
  end 
 if(config2==2)
  Tw_delPwf=input('Enter the value of Tw (in s) for the wash-out circuit [1 - 20]s :  ' );
  TFPSS=KdelPw_pssf.*(((j*w).*Tw_delPwf)./(1+(j*w).*Tw_delPwf)).*((1+(j*w).*T1_delPwf)./(1+(j*w).*T2_delPwf)).*((1+(j*w).*T3_delPwf)./(1+(j*w).*T4_delPwf));
end 
if (config2==3)
    Ks3_delPw=1;
    Tw1_delPwf=input('Enter the value of Tw1 (in s) for the wash-out circuit-1 [1 - 20]s :  ' );
    Tw2_delPwf=input('Enter the value of Tw2 (in s) for the wash-out circuit-2 [1 - 20]s :  ' );
    Tw3_delPwf=input('Enter the value of Tw3 (in s) for the wash-out circuit-3 [1 - 20]s :  ' );
    Tw4_delPwf=input('Enter the value of Tw4 (in s) for the wash-out circuit-4 [1 - 20]s :  ' );
    T6_delPwf=input('Enter the value of T6 (in s) for the slip-path delay [0.01 - 0.05]s :  ' );
    T7_delPwf=input('Enter the value of T7 (in s) for the equivalent integrator [1 - 10]s :  ' );
    T8_delPwf=input('Enter the value of T8 (in s) for the Filter circuit [0 - 0.01]s :  ' );
    T9_delPwf=input('Enter the value of T9 (in s) for the Filter circuit [0.1 - 0.2]s :  ' );
   TFPSS1=(((j*w).*Tw1_delPwf)./(1+(j*w).*Tw1_delPwf)).*(((j*w).*Tw2_delPwf)./(1+(j*w).*Tw2_delPwf)).*(1./(1+(j*w).*T6_delPwf)).*((1+j*w.*T8_delPwf)./(1+(j*w).*T9_delPwf).^2).^4;
   TFPSS2=(((j*w).*Tw3_delPwf)./(1+(j*w).*Tw3_delPwf)).*(((j*w).*Tw4_delPwf)./(1+(j*w).*Tw4_delPwf)).*(((-j*w).*T7_delPwf)./(1+(j*w).*T7_delPwf)).*((((1+j*w.*T8_delPwf)./(1+(j*w).*T9_delPwf).^2).^4).*Ks3_delPw-1);
   TFPSS3=KdelPw_pssf.*((1+(j*w).*T1_delPwf)./(1+(j*w).*T2_delPwf)).*((1+(j*w).*T3_delPwf)./(1+(j*w).*T4_delPwf));
   TFPSS=(TFPSS1+ TFPSS2).*TFPSS3;
end
PHASE2=angle(TFPSS)*180/pi;
PHASEC=PHASE1+PHASE2;
figure(2)
plot(w/2/pi, PHASEC, 'r', w/2/pi, PHASE1, 'b', w/2/pi, PHASE2, 'k');
grid on
clear PHASEC  PHASE2;
mpssd=input('Enter 0 if you are satisfied with the design, otherwise 1 to re-design the pss:  ');
%-------------------------------------------------------------------

end %while loop end if you are satisfied...psstype==2


fprintf('\n Update the delPW_pss.dat for the machine %d \n',rr);
disp('----------------------------------------------------------------------')
disp('gen.no    Ks        T1             T2               T3           T4 ')
disp('-----------------------------------------------------------------------')
fprintf('%d \t %d \t %5.5f \t %5.5f \t  %5.5f \t %5.5f \n\n', rr, KdelPw_pssf, T1_delPwf, T2_delPwf, T3_delPwf, T4_delPwf)
disp('------------------------------------------------------------------------')
end %if(psstype==2) loop end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if(psstype==3)
mpssd=1;
while(mpssd==1)
 delete(get(0,'children'))
    figure(1)
    plot( w/2/pi, PHASE1);
    disp('Enter 1 : Plain power type PSS')
    disp('      2 : With measurement delay and washout time constant')
    config_power = input('Enter your choice : ');
    
    if(config_power==1)
       Kpower_pssf=input('Enter the PSS gain Ks:  ' );
       PHASEC=PHASE1+90;
       figure(2)
       plot(w/2/pi, PHASEC, 'r', w/2/pi, PHASE1, 'b');
       Tw_powerf=10;Td_powerf=0.05;
    end
    
 if(config_power==2)
  Kpower_pssf=input('Enter the PSS gain Ks:  ' );
  Tw_powerf=input('Enter the value of Tw (in s) for the wash-out circuit [1 - 20]s :  ' );
  Td_powerf=input('Enter the value of TR (in s) for the measurement delay [0.01 - 0.05]s :  ' );
  TFPOWER= Kpower_pssf.*(((j*w).*Tw_powerf)./(1+(j*w).*Tw_powerf)).*(1./(1+(j*w).*Td_powerf));
  PHASE2=angle(TFPOWER)*180/pi+90; 
  PHASEC=PHASE1+PHASE2;
  figure(2)
  plot(w/2/pi, PHASEC, 'r', w/2/pi, PHASE1, 'b', w/2/pi, PHASE2, 'k');
 end

grid on
clear PHASEC;
mpssd=input('Enter 0 if you are satisfied with the design, otherwise 1 to re-design the pss:  ');
%-------------------------------------------------------------------

end %while loop end if you are satisfied...psstype==3
fprintf('\n Update the power_pss.dat for the machine %d \n',rr);
disp('----------------------------------------------------------------------')
disp('gen.no     Ks             TR            Tw    ')
disp('-----------------------------------------------------------------------')
fprintf('%d \t %5.5f \t %5.5f \t %5.5f  \n\n', rr, Kpower_pssf, Td_powerf,Tw_powerf)
disp('------------------------------------------------------------------------')
end % if end of psstype==3

disp('Press any key to obtain the amplitude response of GEPS(iw)')
pause
figure(3)
plot(w/2/pi,abs(Te_Vs))
   
clear PHASE1;
clear Xwd;
clear Te_Vs;
    
 




