clear;
clc;
eps=1e-5;
N=0;
A=[25,-41,10,-6;
   -41,68,-17,10;
   10,-17,5,-3;
   -6,10,-3,2];
[m,n]=size(A);
while (1)
   N=N+1;
   B=abs(A);
   max=B(1,2);
   t=1;k=2;
   for i=1:m-1
       for j=i+1:n
           if (B(i,j)>max)
               max=B(i,j);
               t=i;
               k=j;
           end
       end
   end
   if max<eps
       break;
   end
   max=A(t,k);
   fprintf('max=%f,i=%d,j=%d\n',max,t,k);
   if A(t,t)~=A(k,k)
       jz=1/2*(atan(2*A(t,k)/(A(t,t)-A(k,k))));
   elseif A(t,t)==A(k,k)&A(t,k)>0
       jz=pi/4;
   elseif A(t,t)==A(k,k)&A(t,k)<0
       jz=-pi/4;
   end
   fprintf('第%d次旋轉(zhuǎn)角度為：%f\n',N,jz*180/pi);
   A1=A;
   for l=1:n
       A1(t,l)=A(t,l)*cos(jz)+A(k,l)*sin(jz);
   end
   for l=1:n
       A1(l,t)=A(t,l)*cos(jz)+A(k,l)*sin(jz);
   end
   for l=1:n
       A1(k,l)=-A(t,l)*sin(jz)+A(k,l)*cos(jz);
   end
   for l=1:n
       A1(l,k)=-A(t,l)*sin(jz)+A(k,l)*cos(jz);
   end
   A1(t,t)=A(t,t)*(cos(jz)^2)+A(k,k)*(sin(jz)^2)+A(t,k)*sin(2*jz);
   A1(k,k)=A(t,t)*(sin(jz)^2)+A(k,k)*(cos(jz)^2)-A(t,k)*sin(2*jz);
   A1(t,k)=0;
   A1(k,t)=0;
   A=A1;
end
for i=1:n
    Y(i)=A(i,i);
end 
%輸出特征值
Y'