#include ../data/case14.txt
#include ../lib/idx_gen.txt
#include ../lib/idx_bus.txt
#include ../lib/idx_brch.txt
#include ../lib/make_y_bus.txt
#include ../lib/dsbus_dv.txt
#include ../lib/make_sdzip.txt
#include ../lib/make_sbus.txt
#include ../lib/newtonpf.txt
#include ../lib/runpf.txt
#include ../lib/pfsoln.txt
#include ../lib/total_load.txt

info(`Running power flow on case14...`);
V = runpf();
info(`Run power flow on case14 end`);

bus = assign(bus, abs(V), [0], [VM-1, VM] );
bus = assign(bus, angle(V) * 180 / pi, [0], [VA-1, VA] );

on_temp = find(slice(gen, [0], [GEN_STATUS-1, GEN_STATUS]) > 0);
index_gbus_type = get_multi(slice(gen, [0], [GEN_BUS-1, GEN_BUS]), on_temp) - 1;
on = find(get_multi(slice(bus, [0], [BUS_TYPE-1, BUS_TYPE]), index_gbus_type) != PQ)';
off = find(slice(gen, [0], [GEN_STATUS-1, GEN_STATUS]) <= 0)';
gbus = get_multi(slice(gen, [0], [GEN_BUS-1, GEN_BUS]), on)';
Pd_Qd = total_load(select(bus, gbus-1), [], 1, [], [], []);

return Pd_Qd;