%~~~~~~~~~~~~~~~~~~~~~eig_test_pss_final~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%format short g
YES2=input('Enter 1 to display  ALL eigenvalues (EIG), otherwise 0 (using EIGS):  ');
if(YES2)
   [V1 D]=eig(full(AT));
   V1=sparse(V1);
   W=V1\eye(size(V1));
   PF=V1.*conj(W');
   D=diag(D);
   n_eig=1;
   for k=1:size(AT,1)-3
    if (sum(abs(D(k+[0:3])))>0.1)
     n_eig = n_eig+1;
    end
   end
   D=diag(D(1:n_eig))-0.000001;
else
   n_eig=10;
   f_eig=1;
   fprintf('You  are scanning %d eigenvalues around %5.3f Hz .....\n', n_eig,f_eig);
   disp('Please press a key: ')
   pause
   options.tol=1e-12;
   options.disp =0;
   options.maxit = 25;
   [V1 D Flag1]= eigs(AT,n_eig,j*(f_eig*2*pi),options);
   [W1 D Flag2]= eigs((AT)',n_eig,j*(f_eig*2*pi),options);
   VW=conj(V1)'*W1;
   PF=(V1.*W1)/VW;
   EIG=diag(D)-0.000001;
   D=diag(EIG);
   if ((Flag1==1)|(Flag2==1))
      disp('----------------------------NOTE---------------------------------------');
      disp('Determination of eigenvalues with EIGS is NOT accurate. Please use EIG function');
      end
end

SL_lamb=(1:n_eig)';
disp('-------------------------------------------------------------------');
disp(   'SL_number         Eigenvalue       dampingfactor     frequency(Hz)')
disp('--------------------------------------------------------------------');
[SL_lamb diag(D) -real(diag(D))./sqrt(real(diag(D)).*real(diag(D))+imag(diag(D)).*imag(diag(D)))  abs(imag(diag(D)))./(2*pi) ]   
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
disp('---------------------------------------------------------------------');
EIG=diag(D);
rpt=1;
while (rpt==1)
SL_num=input('Enter the serial number of the eigenvalue for which you want to obtain the P.factor:  ');                                                                                              
PF_lamb=PF(:,SL_num);

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

PF_lamb_nr = PF_lamb/max(PF_lamb);
[m_sort m_ind]=sort(abs(PF_lamb_nr));
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (npss~=0)
  if (Tmd_slip_nt==0)
     nm_st_pss_nt = 3;  % if you are using input measuring time delay 
  else
      nm_st_pss_nt = 2; 
  end
  
  if (Tmd_slip_t==0)
     nm_st_pss_t = 5;  % if you are using input measuring time delay 
  else
      nm_st_pss_t = 4; 
   end
end %-----End of npss if loop-------------------


 rev_ind=[size(AT,1):-1:1];

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
n_ind=m_ind(rev_ind);
n_ind=n_ind+1;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
X1=('');
X3=('');
for i=1:ngen
X1=strvcat(X1,'Delta-','Slip-','Field-','DampH-','DampG-','DampK-','Efd-','VR_DC-','xB_DC_AC-','xF_DC-','x1_st_tu-','x2_st_tu-','x3_st_tu-','y1_gv-','PG-','z-');
X2d=gen(i,1)*ones(16,1);
X2=num2str(X2d);
X3=strvcat(X3,X2);
end

Xt=strcat(X1,X3);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if(npss~=0)
   XPt_nt =('');
   XPt_t=('');
   XP_nt=('');
   X4_nt=('');
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if(sum(xpss_nt)~=0)
for i=1:size(xpss_nt,2);
    
  if (Tmd_slip_nt==0)
     XP_nt=strvcat(XP_nt,'xpss_1-','xpss_2-','xpss_3-');  % if you are using input measuring time delay 
  else
      XP_nt=strvcat(XP_nt,'xpss_1-','xpss_2-'); 
  end

   XPd_nt=xpss_nt(i)*ones(nm_st_pss_nt,1);
   XP2_nt=num2str(XPd_nt);
   X4_nt=strvcat(X4_nt,XP2_nt);
end
XPt_nt=strcat(XP_nt,X4_nt);
end
     %-------------------------------------------------
if(sum(xpss_t)~=0)
   XP_t=('');
   X4_t=(''); 
 for i=1:size(xpss_t,2);
     
  if (Tmd_slip_t==0)
     XP_t=strvcat(XP_t,'xpss_1-','xpss_2-','xpss_3-','xpss_4-','xpss_5-');  % if you are using input measuring time delay 
   else
      XP_t=strvcat(XP_t,'xpss_1-','xpss_2-','xpss_3-','xpss_4-');
   end
  
   XPd_t=xpss_t(i)*ones(nm_st_pss_t,1);
   XP2_t=num2str(XPd_t);
   X4_t=strvcat(X4_t,XP2_t);
  end
XPt_t=strcat(XP_t,X4_t);
end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 XPt=strvcat(XPt_nt,XPt_t);
 Xt=strvcat(Xt,XPt);
end   %-----End of npss if loop-------
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if(npss1~=0)
   XdelPW=('');
   X5=('');
for i=1:size(delPw_pss,1);
   XdelPW=strvcat(XdelPW,'xdelPw_1-','xdelPw_2-','xdelPw_3-','xdelPw_4-','xdelPw_5-','xdelPw_6-','xdelPw_7-','xdelPw_8-','xdelPw_9-','xdelPw_10-','xdelPw_11-','xdelPw_12-','xdelPw_13-','xdelPw_14-','xdelPw_15-','xdelPw_16-');
   XdelPWd=delPw_pss(i,1)*ones(16,1);
   XdelPWd2=num2str(XdelPWd);
   X5=strvcat(X5,XdelPWd2);
end
  XPt=strcat(XdelPW,X5);
  Xt=strvcat(Xt,XPt);
end  %-----End of npss1 if loop-------

if(npss2~=0)
   Xpower=('');
   X6=('');
for i=1:size(power_pss,1);
   Xpower=strvcat(Xpower,'xpower_1-','xpower_2-');
   Xpowerd=power_pss(i,1)*ones(2,1);
   Xpowerd2=num2str(Xpowerd);
   X6=strvcat(X6,Xpowerd2);
end
  XPt=strcat(Xpower,X6);
  Xt=strvcat(Xt,XPt);
end  %-----End of npss2 if loop-------
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Yt=Xt(n_ind,:);
PF_lamb_sort=PF_lamb_nr(n_ind-1);

disp('----------------------------------------------------------------------------------');
disp('      State variable      Mag(Norm PF)    ang(Norm PF)deg.  Mag(PF)    ang(PF)deg.');
disp('----------------------------------------------------------------------------------');
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   for i=1:(size(AT,1))
      if (abs(PF_lamb_sort(i))>=0.1)
      fprintf('\t %s %15.4f %13.2f %15.4f %13.2f\n',Yt(i,:),abs(PF_lamb_sort(i)),angle(PF_lamb_sort(i))*180/pi,abs(PF_lamb_sort(i)*full(max(PF_lamb))),angle(PF_lamb_sort(i)*full(max(PF_lamb)))*180/pi);
       end
    end


%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 nonosc=EIG(SL_num);
if(abs(imag(nonosc))<1e-3)
      disp('-----------------------------------------------------------------------------');
      disp('It is a NON-OSCILLATORY MODE...')
      disp('-----------------------------------------------------------------------------');
else
 slip_PF_lamb=PF_lamb(1:16:ngen*16);
 slip_PF_nr=slip_PF_lamb/max(slip_PF_lamb);
[m_sort_slip  m_ind_slip]=sort(abs(slip_PF_nr));
rev_ind_slip=[ngen:-1:1];
n_ind_slip=m_ind_slip(rev_ind_slip);
n_sort_slip=m_sort_slip(rev_ind_slip);
gen_num=gen(n_ind_slip,1);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
abs_sort=0.05;
k=1;
r_sort_slip(k)=0;
for i=1:ngen
   if(abs(n_sort_slip(i))>abs_sort)
      r_sort_slip(k)=n_sort_slip(i);
      k=k+1;
   end
end
r_ind_slip=n_ind_slip(1:size(r_sort_slip,2));
par_gen_numd=gen_num(1:size(r_sort_slip,2));
   
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
V=V1(:,SL_num);
slip_V=V(1:16:ngen*16);
%-----To filter out a mode in which participation of slip is very small----
splip_PF_norm = sparse(zeros(1,nb));
splip_PF_norm(gen(:,1)) = PF_lamb_nr(1:16:ngen*16);
splip_PF_norm_filt = splip_PF_norm(par_gen_numd);
slip_Vpart = 0.3;     %  use this variable to decide the slip part. level
if(sum(abs(splip_PF_norm_filt))< slip_Vpart)
    disp('---------------------------------------------------------------------------------');
    disp('You have chosen a NON-SWING MODE with very low slip participation...')
    disp('---------------------------------------------------------------------------------');
  else
n_slip_Vd=slip_V(r_ind_slip);
n_slip_Vd_nr=n_slip_Vd/max(n_slip_Vd);
[m_sort_slipV  m_ind_slipV]=sort(abs(n_slip_Vd_nr));
rev_ind_slipV=[size(n_slip_Vd,1):-1:1];
n_ind_slipV=m_ind_slipV(rev_ind_slipV);
par_gen_num=par_gen_numd(n_ind_slipV);
n_slip_V=n_slip_Vd(n_ind_slipV);
slip_V_angle=angle(n_slip_V)*180/pi;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
swing_ang=90;
p1=1;
p2=1;
Group1(p1)=par_gen_num(1);
Group2(p2)=0;
a=slip_V_angle(1);
for i=2:size(r_sort_slip,2)
   b=slip_V_angle(i);
   if((abs(a-b)>swing_ang)&(abs(a-b)<(360-swing_ang )))
      Group2(p2)=par_gen_num(i);
      p2=p2+1;
   else
      p1=p1+1;
      Group1(p1)=par_gen_num(i);
   end
end
%------------------------------------------------------------------
disp('----------------------------------------------------------------------------------');
if(sum(Group2)~=0)
   disp('      You have chosen a SWING-MODE       ')
else
   disp('      You have chosen a NON SWING-MODE   ')
end
%-------------------------------------------------------------------


disp('----------------------------------------------------------------------------------');
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

disp('The generator(s) in group-1 is(are) ...       ')
Group1
disp('The generator(s) in group-2 is(are) ...      ')
Group2
disp('----------------------------------------------------------------------------------');

%--------------------------------------------------------------------------

TRACE=input('Enter 1 if you want to plot the compass plot, otherwise 0:   ');
if(TRACE)
  disp('NOTE: Use mouse click on the plot to identify the generator')
  disp('Press any key')
  pause
  r_eig_plot
end
end %It is non swing mode of less participation if loop end
end %It is Non-oscillatory mode if loop end.
rpt=input('Enter 1 if you want to repeat  for another eigenvalue, otherwise 0: ') ;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
clear Group1 Group2;
clear slip_V_angle;
clear n_slip_Vd;
clear n_slip_V;
clear par_gen_numd;
clear par_gen_num;
clear r_ind_slip;
clear rev_ind;
clear n_ind;
clear p1 p2;
clear r_sort_slip k;
clear gen_num;
clear rev_ind_slipV;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
end

