clc
clear all 

%% params
wn = 1;
zn = 0.005;

frf = @(w) 1./((wn^2-w.^2)+1j*(2*zn*wn*w));

%% Data
Nw = 100;
Ws = linspace(0, 2, Nw)';
Hs = frf(Ws);

%% Loewner Matrices
is = [20 40];

mus = Ws(is(1))*[1j; -1j];
v = [Hs(is(1)); conj(Hs(is(1)))];

lams = Ws(is(2))*[1j; -1j];
w = [Hs(is(2)); conj(Hs(is(2)))];

R = ones(1, 2);
C = ones(2, 1);

L = (v*R-C*transpose(w))./(mus*R-C*transpose(lams));
LS = ((mus.*v)*R-C*transpose(lams.*w))./(mus*R-C*transpose(lams));

Jmu = 1/sqrt(2)*[1 -j;1 j];
Jlam = Jmu;

sys = ss(-LS, v, transpose(w), 0);
sys.E = -L;

sysj = ss(-real(Jmu'*LS*Jlam), real(Jmu'*v), real(transpose(w)*Jlam), 0);
sysj.E = -real(Jmu'*L*Jlam);

[Hmk, Hak] = bode(sys, Ws); 
[Hmjk, Hajk] = bode(sysj, Ws);


figure(1);
plot(Ws, squeeze(Hmk), '.-', Ws, squeeze(Hmjk), 'o-', Ws, abs(squeeze(Hs)), 'k.-', ...
  abs([mus; lams]), abs([v; w]), 'k*')