% eigen2.m  demonstrate more eigenvalue testing
clear
format compact
sprintf('output file eigen2_m.out being written \n')
fid = fopen('eigen2_m.out', 'w');
fprintf(fid, 'eigen2.m running\n'); 
e1=1+j*2; % desired eigenvalues, complex
e2=2+j*3;
e3=3+j*4;
e4=4+j*5;
fprintf(fid,'desired eigenvalues= \n');
fprintf(fid,'%f+j%f  %f+j%f %f+j%f %f+j%f \n', real(e1),imag(e1), real(e2),imag(e2), real(e3), imag(e3), real(e4),imag(e4));
P=poly([e1 e2 e3 e4]) % make polynomial from roots
fprintf(fid,'\npolynomial coefficients= \n');
fprintf(fid,'%f+j%f  %f+j%f %f+j%f \n%f+j%f %f+j%f\n', real(P(1)),imag(P(1)), real(P(2)),imag(P(2)), real(P(3)), imag(P(3)), real(P(4)),imag(P(4)), real(P(5)),imag(P(5)));

r=roots(P) % check that roots come back
fprintf(fid,'\nroots of polynomial = desired eigenvalues= \n');
fprintf(fid,'%f+j%f  %f+j%f %f+j%f %f+j%f \n',real(r(1)),imag(r(1)),real(r(2)),imag(r(2)),real(r(3)),imag(r(3)),real(r(4)), imag(r(4)));

A=zeros(4);
A(1,1)=-P(2); % build matrix A
A(1,2)=-P(3);
A(1,3)=-P(4);
A(1,4)=-P(5);
A(2,1)=1.0;
A(3,2)=1.0;
A(4,3)=1.0;
fprintf(fid,'\nA matrix = \n');
for i=1:4
  fprintf(fid,'%f+j%f  %f+j%f %f+j%f %f+j%f \n', real(A(i,1)),imag(A(i,1)), real(A(i,2)),imag(A(i,2)), real(A(i,3)), imag(A(i,3)), real(A(i,4)),imag(A(i,4)));
end

[v e]=eig(A) % compute eigenvectors and eigenvalues
fprintf(fid,'\ncomputed eigen values= \n');
fprintf(fid,'%f+j%f  %f+j%f %f+j%f %f+j%f \n', real(e(1,1)),imag(e(1,1)), real(e(2,2)),imag(e(2,2)), real(e(3,3)), imag(e(3,3)), real(e(4,4)),imag(e(4,4)));
fprintf(fid,'\ncomputed eigen vectors = \n');
for i=1:4
  fprintf(fid,'%f+j%f  %f+j%f %f+j%f %f+j%f \n', real(v(i,1)),imag(v(i,1)), real(v(i,2)),imag(v(i,2)), real(v(i,3)), imag(v(i,3)), real(v(i,4)),imag(v(i,4)));
end
I=eye(4);    % identity matrix
maxerr = 0.0;
for i=1:4
  err=det(A-I.*e(i,i)); % should be near zero
  maxerr = max(maxerr,abs(err));
end
fprintf(fid,'\ndet(A-Ie) maxerr= %f \n',maxerr);
maxerr = 0.0;
for i=1:4
  err=A*v(:,i)-e(i,i)*v(:,i); % note columns, near zero
  maxerr = max(maxerr,abs(err));
end
fprintf(fid,'\nA*v(:,i)-e(i,i)*v(:,i) maxerr= %f \n',maxerr);
z3=trace(A)-trace(e); % should be near zero
fprintf(fid,'\ntrace(A-trace(r)= %f \n',z3);

fprintf(fid,'\ncompute B=Q*A*inv(Q) and eigenvalues of B should be same \n');
Q=randn(4);
B=Q*A*inv(Q);
B=B.*(5+j*6); % scale then unscale eigenvalues
[v e]=eig(B); % compute eigenvectors and eigenvalues
ee=e./(5+j*6); % expect e1, e2, e3, e4
fprintf(fid,'computed eigenvalues= \n');
fprintf(fid,'%f+j%f  %f+j%f %f+j%f %f+j%f \n', real(ee(1,1)),imag(ee(1,1)), real(ee(2,2)),imag(ee(2,2)), real(ee(3,3)), imag(ee(3,3)), real(ee(4,4)),imag(ee(4,4)));

I=eye(4);    % identity matrix
maxerr = 0.0;
for i=1:4
  err=det(B-I.*e(i,i)); % should be near zero
  maxerr = max(maxerr,abs(err));
end
fprintf(fid,'\ndet(B-Ie) maxerr= %f \n',maxerr);
maxerr = 0.0;
for i=1:4
  err=B*v(:,i)-e(i,i)*v(:,i); % note columns, near zero
  maxerr = max(maxerr,abs(err));
end
fprintf(fid,'\nB*v-e*v maxerr= %f \n',maxerr);
z3=trace(B)-trace(e); % should be near zero
fprintf(fid,'\ntrace(B-trace(r)= %f \n',z3);
