function [r,J] = loadFlow(x,sref,z,p)
%load flow residuals
%[r,J] = loadFlow(x,sref,z,p)
%inputs:
%x: network voltages; x = [V2R V2I V3R V3I]; 
%sref: network injections; sref = [P2+1i*Q2 P3+1i*Q3];
%z: dynamic parameters
%p: static parameters
%outputs
%r: residual; 
%J: Jacobian of the residual; J = dr/dx

X1 = p(1);
X2 = p(2);
n = z(1);

%slack bus voltage is used for reference
V0R = 1;
V0I = 0;
%retrieve the remaining voltage guesses
V2R = x(1);
V2I = x(2);
V3R = x(3); 
V3I = x(4); 
%retrieve the actual injections
S2r = sref(1);
S3r = sref(2);

%reconstruct complex states to simplify residue expressions
V0 = V0R + V0I*1i;
V2 = V2R + V2I*1i;
V3 = V3R + V3I*1i;

%impedences
Z1 = 1i*X1;
Z2 = 1i*X2;

%compute the injections from the voltage guesses
S2 = n/conj(Z1)*conj(V0-n*V2)*V2 - conj(V2-V3)/conj(Z2)*V2;
S3 = conj(V2-V3)/conj(Z2)*V3;

%compute the residuals
r = zeros(4,1);
r(1) = real(S2- S2r);
r(2) = imag(S2- S2r);
r(3) = real(S3- S3r);
r(4) = imag(S3- S3r);

%compute the Jacobian
dS2dV2 = n/conj(Z1)*(conj(V0 - n*V2)*[1 1i] - n*[1 -1i]*V2) - 1/conj(Z2)*(conj(V2-V3)*[1 1i] + [1 -1i]*V2);
dS2dV3 = 1/conj(Z2)*[1 -1i]*V2;


dS3dV2 = 1/conj(Z2)*[1 -1i]*V3;
dS3dV3 = 1/conj(Z2)*(-[1 -1i]*V3 + conj(V2-V3)*[1 1i]);

J = [...
    real(dS2dV2) real(dS2dV3);
    imag(dS2dV2) imag(dS2dV3);
    real(dS3dV2) real(dS3dV3);
    imag(dS3dV2) imag(dS3dV3);
    ];


IspIc = conj(S2r/V2);






