function [Hf, Ht] = makePTDFac(baseMVA, bus, branch, Pslack, Vfixed, Qslack)
%MAKEPTDFAC   Builds the AC PTDF matrix.
%   H = MAKEPTDFAC(BASEMVA, BUS, BRANCH, PSLACK, VFIXED, QSLACK)
%   returns the AC PTDF matrix. The matrix is 2*nbr x 2*nb
%   where nbr is the number of branches and nb is the number of buses.
%   PSLACK is the bus index of the slack bus for real power, VFIXED is a
%   vector of bus indices where the voltage magnitude is assumed to remain
%   fixed, and QSLACK is a vector of indices of buses whose Q injections are
%   assumed to be a function of the other injections (uses VFIXED as default).
%
%   Examples:
%       H = makePTDFac(baseMVA, bus, branch);
%       H = makePTDFac(baseMVA, bus, branch, Pslack);
%       H = makePTDFac(baseMVA, bus, branch, Pslack, Vfixed);
%       H = makePTDFac(baseMVA, bus, branch, Pslack, Vfixed, Qslack);

%   MATPOWER
%   Copyright (c) 2006-2015 by Power System Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See http://www.pserc.cornell.edu/matpower/ for more info.

[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;

if nargin < 6
    if nargin < 5
        Vfixed = find(bus(:, BUS_TYPE) == PV | bus(:, BUS_TYPE) == REF);
        if nargin < 4
            Pslack = find(bus(:, BUS_TYPE) == REF);
        end
    end
    Qslack = Vfixed;
end

nb = size(bus, 1);
nbr = size(branch, 1);
allbuses = (1:nb)';
ref     = find(bus(:, BUS_TYPE) == REF); ref = ref(1);
noPslack = find(allbuses ~= Pslack);
noQslack = allbuses; noQslack(Qslack) = 0; noQslack = find(noQslack);
noref   = find(allbuses ~= ref);
noVfixed = allbuses; noVfixed(Vfixed) = 0; noVfixed = find(noVfixed);
noslack = [ noPslack; nb+noQslack ];
norefVfixed = [ noref; nb+noVfixed ];

%% check that bus numbers are equal to indices to bus (one set of bus numbers)
if any(bus(:, BUS_I) ~= (1:nb)')
    error('makePTDFac: buses must be numbered consecutively in bus matrix; use ext2int() to convert to internal ordering')
end
if length(Vfixed) ~= length(Qslack)
    error('makePTDFac: Vfixed and Qslack must have the same dimensions')
end

%% compute PTDF
[Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
V = bus(:, VM) .* exp(sqrt(-1) * bus(:, VA) * pi/180);
[dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V);
[dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St] = dSbr_dV(branch, Yf, Yt, V);
Bbus = [    real(dSbus_dVa) real(dSbus_dVm);
            imag(dSbus_dVa) imag(dSbus_dVm) ];
Bf = [  real(dSf_dVa) real(dSf_dVm);
        imag(dSf_dVa) imag(dSf_dVm) ];
Bt = [  real(dSt_dVa) real(dSt_dVm);
        imag(dSt_dVa) imag(dSt_dVm) ];
Hf = zeros(2*nbr, 2*nb);
Ht = zeros(2*nbr, 2*nb);
Bbusinv = inv(Bbus(noslack, norefVfixed));
Hf(:, noslack) = full(Bf(:, norefVfixed) * Bbusinv);
Ht(:, noslack) = full(Bt(:, norefVfixed) * Bbusinv);
