function [JxSystem,JdSystem] = jfun(t,x,d,z,p,St,Co,fnames)
%system Jacobian function
%z = vector of discrete states
%p = vector of parameter values
%St = assignment of states to modules
%Co = connectivity matrix defnining shared states
%modules: 
%   - module 1: core network
%   - module 2: timer transformer
%   - module 3: timer capacitor
%   - module 4: load
%   - module 5: capacitor bank
%   - module 6: wind turbine

nMo = size(St,1);
nSt = St(end);
JxSystem = zeros(nSt);
JdSystem = zeros(nSt);
ni = 0;

for i = 1:nMo
    xCur = x(St(i,1):St(i,2));
    dCur = d(St(i,1):St(i,2));
    [JxCur,JdCur] = fnames{i,3}(t,xCur,dCur,z{i},p{i});
    nCur = size(JxCur,1);
    JxSystem(ni+1:ni+nCur,St(i,1):St(i,2)) = JxCur;
    JdSystem(ni+1:ni+nCur,St(i,1):St(i,2)) = JdCur;
    ni = ni + nCur;
end
c = size(Co,1);
for i = 1:c
    ind1 = St(Co(i,1),1) - 1 + Co(i,2);
    ind2 = St(Co(i,3),1) - 1 + Co(i,4);
    JxSystem(ni+i,ind1) = 1;
    JxSystem(ni+i,ind2) = -1;
end
%states with explicit derivatives: V0-(1,2), xTT-(13), xTC-(15), xd-(17), wr-(30) 
%residues of explicit derivatives: V0-(5,6), xTT-(7), xTC-(9), xd-(11), wr-(20)
% vecxx = [3:12,14,16,18:29,31,32];
% vecxr = [1:4,8,10,12:19,21:32];
% vecdx = [1,2,13,15,17,30];
% vecdr = [5,6,7,9,11,20];
% svd(JxSystem(vecxr,vecxx))
% eig(JdSystem(vecdr,vecdx))
% % svd(JdSystem)
% pause