% SEA-COND方法（SEA為 Sensor Elimination Algorithm, 即傳感器淘汰算法，COND為條件數, Condition number）：
 % 直接計算模態矩陣 A 的條件數（Condition number）, cond(A)
 % 模態矩陣 A 相當于是由候選傳感器“測量”的模態矩陣。
 % 當前要淘汰的傳感器（自由度）是按下列方式確定的：
 % 輪流刪除每個傳感器，即模態矩陣 A 中的每一行，然后，計算模態矩陣 A 的條件數 cond(A)，
 % 哪一種刪除所產生的條件數cond(A)值最小，那么，相應的傳感器就是本輪要淘汰的對象。
 
clear;
%============================ can be updated by user===================================================================
modal=[0.0157    0.0302   -0.0422   -0.0506    0.0545    0.0532   -0.0463    0.0343    0.0183
       0.0299    0.0488   -0.0496   -0.0313   -0.0000   -0.0329    0.0545   -0.0556   -0.0348
       0.0411    0.0488   -0.0161    0.0313   -0.0545   -0.0329   -0.0177    0.0556    0.0479
       0.0483    0.0302    0.0306    0.0506    0.0000    0.0532   -0.0337   -0.0343   -0.0563
       0.0508   -0.0000    0.0521   -0.0000    0.0545    0.0000    0.0573    0.0000    0.0592
       0.0483   -0.0302    0.0306   -0.0506    0.0000   -0.0532   -0.0337    0.0343   -0.0563
       0.0411   -0.0488   -0.0161   -0.0313   -0.0545    0.0329   -0.0177   -0.0556    0.0479
       0.0299   -0.0488   -0.0496    0.0313   -0.0000    0.0329    0.0545    0.0556   -0.0348
       0.0157   -0.0302   -0.0422    0.0506    0.0545   -0.0532   -0.0463   -0.0343    0.0183];
      %Modal Matrix中振型向量按列放置，不要求一定是方陣
n=9;          %n-measured DOF,
m1=2;m2=2;    %m1,m2-modal order numbers between which the mode shapes are concerned, set by user
%===================================================================================================================
m=m2-m1+1;    %m-number of concerned modes 
a(1:n,1:m)=modal(1:n,m1:m2);
nj=n;
priority(1:n)=0;

for j=1:1:n %輪次
   for i=1:1:nj %輪流刪除每個自由度,同時計算每種刪除所產生的行列式的值.刪掉第i自由度的方式是將第i行從a中擠出，新矩陣記為ai
      for k=1:1:i-1
         ai(k,1:m)=a(k,1:m);    %i行以前不變
      end
      for k=i:1:nj-1             %i行以后均提前一行
         ai(k,1:m)=a(k+1,1:m);
      end
      d(i)=cond(ai,2);         %計算模態矩陣 Ai 的條件數
   end
   d  %本輪完成后，顯示各淘汰方式的條件數的值
   %找出本輪中各條件數的最小值，然后記下其對應的自由度號
   dmin=1e100;imin=0;
   for i=1:1:nj
      if d(i)<dmin
         dmin=d(i);
         imin=i;
      end
   end
   %將本輪結果填入設置傳感器的優先級向量priority
   k=0;
   for i=1:1:n
      if priority(i)==0
         k=k+1;
         if k==imin
            priority(i)=j/n;   %第i個自由度在第j輪被淘汰，于是其被設置傳感器的優先級被定義為j/n，此值越大，優先布置傳感器的級別越高
            dele(j)=i;      %本輪刪除第i個自由度，dele 僅為用于數字顯示淘汰順序，可以不要此行
         end
      end
   end
   
   %根據上述最小值，確定下一輪的模態矩陣
   for k=1:1:imin-1
       aj(k,1:m)=a(k,1:m);
   end
   for k=imin:1:nj-1
       aj(k,1:m)=a(k+1,1:m);
   end
   a=aj;
   nj=nj-1;
end
figure(1)
bar(priority(1:n))%傳感器的優先級
grid on
xlabel('No. of DOF','fontname','times new roman','fontsize',10,'fontweight','bold');
ylabel('Priority of Instrumentation','fontname','times new roman','fontsize',10,'fontweight','bold');
colormap([1 0 0]);
dele%數字顯示淘汰順序
figure(2)%淘汰順序
DeleMatrix(1:n,1:n)=1;%
for i=2:1:n
   for j=1:1:i-1
       DeleMatrix(i,dele(j))=0;
   end
end
bar3(DeleMatrix);
DeleMatrix