
% fea_transmissibility_core_complex.m  ver 1.1 by Tom Irvine

function[f,acc,acc_complex,rd,rd_complex]=fea_transmissibility_core_complex(f,nff,nem,ModeShapes,MST,...
                             Mwd,damp,omegan,TT,ngw,TZ_tracking_array,node,num_modes)

ij   = TZ_tracking_array(node);

np   = numel(f);
omega = 2*pi*f(:);                 % column
om2   = omega.^2;

% modal data (columns)
damp   = damp(:);
omegan = omegan(:);
omn2   = omegan.^2;

% excitation vector for A
y  = ones(nem,1);
A  = -MST*Mwd*y;                   % nff×1

acc = zeros(np,1,'like',omega);
rd  = zeros(np,1,'like',omega);

acc_complex = zeros(np,1,'like',omega);
rd_complex  = zeros(np,1,'like',omega);

% how many modes to use this pass
nffx = min([nff, numel(damp), numel(omegan), numel(om2), num_modes]);
idx  = 1:nffx;

% precompute column-shaped pieces
Ai    = A(idx);
omn2i = omn2(idx);
itdi  = (2*1i) * (damp(idx).*omegan(idx));    % nffx×1

% reduce ModeShapes to active columns
Md = ModeShapes(:, idx);

% map desired full-DOF index once
pos = find(ngw == ij, 1);
assert(~isempty(pos), 'Index ij not found in ngw.');

oneNem = ones(nem,1,'like',A);

progressbar;
for k = 1:np
    progressbar(k/np);

    den  = omn2i - om2(k) + omega(k).*itdi;   % nffx×1
    Xk   = Md * (Ai ./ den);                  % equals ModeShapes*N with N(idx)=Ai./den

    Ud_k = (-1/om2(k)) * oneNem;

    U    = TT * [Ud_k; Xk];

    uij  = U(pos);


    % absolute acceleration (complex) from displacement phasor uij:
    acc_complex(k) = -om2(k) * (uij);   % <-- add minus sign here
    
    % relative displacement (complex) vs. base (your definition):
    rd_complex(k)  = (uij - Ud_k(1));
    
    % magnitudes if you still want them plotted
    acc(k) = abs(acc_complex(k));
    rd(k)  = abs(rd_complex(k));


end
progressbar(1);
pause(0.2);
