function [Jflow, H, Pf, Qf] = ...
    getLinFlows(V0,A0, n, m, Yf, J_full, mpc)
% ---Jflow is the sensitivity matrix showing line flow change to the change in
% voltage angle and magnitude---
% --H is the AC PTDF matrix---
% ---Pf and Qf are the active and reactive power flow on lines at the current
% operating point---
% ------------ code assumes that only one slack bus (id = 1) is present and no PV bus, 
%so if you wish to generalize the code for the arbitrary number of slack
%and PV buses, please change accordingly Jacobian (J_full and line 34-40)
%of the code to make it more generic. --------------------

define_constants;
% Abus(mpc.gen(:,GEN_BUS)) = mpc.gen(:,VA);
% Sbus = makeSbus(mpc.baseMVA, mpc.bus, mpc.gen);
Cf= getCf(mpc);

%%%% Power Flow Linearization Terms
V = V0.*exp(1j*A0);
Vf = Cf*V;
If = Yf*(V);
E = (diag(1./V0)*V);
Sf = diag(Vf)*conj(If);
Pf = real(Sf); Qf = imag(Sf);
dSf_dA1 = diag(conj(If))*Cf*diag(V);
dSf_dA2 = diag(Vf)*(conj(Yf) * conj(diag(V)));
dSf_dA = 1j*(dSf_dA1 - dSf_dA2);
dSf_dV1 = diag(conj(If))*Cf*diag(E);
dSf_dV2 = diag(Vf)*conj(Yf)*diag(conj(E));
dSf_dV = dSf_dV1 + dSf_dV2;
dP_dAdV = [real(dSf_dA) real(dSf_dV)];
dQ_dAdV = [imag(dSf_dA) imag(dSf_dV)];
Aflow = [dP_dAdV;dQ_dAdV]; 

Jflow = [Aflow(:,2:n) Aflow(:,n+2:n*2)];
J_red = [[J_full(2:n,2:n) J_full(2:n,n+2:n*2)]; [J_full(n+2:n*2,2:n) J_full(n+2:n*2,n+2:n*2)]];
% inv_Jac = inv(J_red);
H = zeros(m*2,n*2);
H_temp = Jflow/J_red;
H(:,2:n) = H_temp(:,1:n-1);
H(:,n+2:n*2) = H_temp(:,n:end);

