% stitest
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

n = 128; nH = 4; nA = 4;
w = 2*pi*linspace(1,5,n)';
p1 = -2 + 20*j; p1cc = -2 - 20*j; m1 = 1;
p2 = -2 + 15*j; p2cc = -2 - 15*j; m2 = 0; 
H = p1 ./ (p1 - j*w) * exp(2*pi*j*m1/nH*(0:nH-1)) + ...
   p1cc ./ (p1cc - j*w) * exp(-2*pi*j*m1/nH*(0:nH-1)) + ...
    p2 ./ (p2 - j*w) * exp(2*pi*j*m2/nH*(0:nH-1)) + ...
   p2cc ./ (p2cc - j*w) * exp(-2*pi*j*m2/nH*(0:nH-1));
%rand('normal')
Hn = H + 2*(rand(n,nH) + j*rand(n,nH));
[A,B,RAB,h] = sti(w,Hn,nA,nA);
plot(w,[real(H) imag(H) real(Hn) imag(Hn)],'-r',w,[real(h) imag(h)],'-w');
pause

[r,p,d,Rrpd] = rpdcov(A,B,RAB);
[x,Rx] = cpaircov([p;r(:)],Rrpd);

axis('square')
axis([-30 30 -30 30])
plot(0,0,'.w'), hold on
pause

cir = exp(j*linspace(0,2*pi,100));
k = 0;
for kA = 1:nA
 for kH = 1:nH
  k = k+1;
  if imag(p(kA)) == 0
   plot(r(kA,kH),0,'+r')
   plot(r(kA,kH)+sqrt(Rx(nA+k,nA+k))*[-1 1],[0 0],'-r')
  elseif imag(p(kA)) >= 0
   [V,D] = eig(Rx(k+nA+(0:1),k+nA+(0:1)));
   temp = [1 j] * (V*(sqrt(D) * [real(cir);imag(cir)]));
   plot(r(kA,kH),'+r')
   plot(r(kA,kH) + temp, '-r')
  end
 end
end
hold off
axis('normal')