 YES2=input('Enter 1 to use EIG function, otherwise 0 to use EIGS function:  ');
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 using EIGS function.\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)')
[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('---------------------------------------------------------');
disp('---------------------------------------------------------');
EIG=diag(D);
nodp=0; % to count number of swing modes with damping factor less than 0.05  
disp('swing modes for which dampingfactor is less than 0.05')  
disp('---------------------------------------------------------'); 
disp(   'SL_number          Eigenvalue    dampingfactor       frequency(Hz)')

for dp=1:1:n_eig-1
   if ((abs(imag(EIG(dp)))>1e-3)&(abs(imag(EIG(dp)))~=abs(imag(EIG(dp+1))))) %To filter out Non_Oscillatory modes 
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
PF_lamb=PF(:,dp);
PF_lamb_nr = PF_lamb/max(PF_lamb);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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;% Sorting Slips for which Normalised slip_PF_nr  is greater than 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(:,dp);
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)
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; %Angle Calculations..
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
damp_fac_lim = 0.05;
if(sum(Group2)~=0) % Only swing modes, Filteringout Non_Swingmodes
     if (-real(EIG(dp))./sqrt(real(EIG(dp)).*real(EIG(dp))+imag(EIG(dp)).*imag(EIG(dp)))<damp_fac_lim) 
         %swing modes to which dp is less than 0.05
      [dp   EIG(dp)   -real(EIG(dp))./sqrt(real(EIG(dp)).*real(EIG(dp))+imag(EIG(dp)).*imag(EIG(dp)))   abs(imag(EIG(dp)))./(2*pi)]
      nodp=nodp+1;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
X1=('');
X3=('');
for i=1:size(r_sort_slip,2)
X1=strvcat(X1,'Slip-');
X2d=par_gen_numd(i,1);
X2=num2str(X2d);
X3=strvcat(X3,X2);
end

Xt=strcat(X1,X3);

slip_PF_sortd=slip_PF_nr(r_ind_slip);

if(npss~=0)
    disp('PSS is Enabled......')
else
    disp('PSS is Disabled.....')
end
disp('---------------------------------------------------------');
disp('State variable     Mag(slip-PF-nr)       angle(slip-PF-nr) in deg.');
disp('---------------------------------------------------------');
 
for i=1:size(r_sort_slip,2)
   fprintf('%s \t\t %12.4f \t\t %10.2f \n',Xt(i,:),abs(slip_PF_sortd(i)),angle(slip_PF_sortd(i))*180/pi);
end 
      %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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;
clear Xt;
clear slip_PF_sortd;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

     end  %swing modes to which dp is less than 0.05 loop end

     end  % (Group2(p2)~=0) if loop end..
   
  end % (sum(abs(splip_PF_norm_filt))> slip_Vpart) if loop end
  
end %(imag(D(dp,dp))>1e-3)
end % for loop end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

disp('---------------------------------------------------------')
if(nodp==0)
    disp('NO swing modes contains damping factor less than 0.05')
else
rpt=1;
while(rpt==1)                                                                                       
disp('Please press a key to obtain the angle of GEPS(s) for the selected machine ')
pause
pss_design %PSS Design...........
rpt=input('Enter 1 if you want to design PSS  for another machine, otherwise 0: ') ;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
end  % while loop end............
  %else 
  %disp('All Eigenvalues contained dampingfactor of greater than 0.05,no need of design pss')
end % if loop end  if nodp loop

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;
clear Xt;
clear slip_PF_sortd;