function varargout = runpfbyopf(casedata, mpopt, fname, solvedcase)
%RUNPFBYOPF  Runs a power flow using runopf.
%   [RESULTS, SUCCESS] = RUNPFBYOPF(CASEDATA, MPOPT, FNAME, SOLVEDCASE)
%% default arguments
if nargin < 4
    solvedcase = '';                %% don't save solved case
    if nargin < 3
        fname = '';                 %% don't print results to a file
        if nargin < 2
            mpopt = mpoption;       %% use default options
            if nargin < 1
                casedata = 'case33bw'; %% default data file is 'case33bw.m' is a radial case
            end
        end
    end
end

define_constants;
mpc = loadcase(casedata);
mpc.branch(33:37,BR_STATUS)=0;       %opening tie lines present
%% now I want to add a generator at bus 10 with reactive power support (without supplying real power)to get constant bus voltage(0.94) at bus 33
mpc.gen(2,:) = [10	0	5	20	-20	1	100	1	10	0	0	0	0	0	0	0	0	0	0	0	0]


mpci = ext2int(mpc);
[bus, gen] = deal(mpci.bus, mpci.gen);
%% save original limits
Vmin = mpc.bus(:, VMIN);
Vmax = mpc.bus(:, VMAX);
Pmin = mpc.gen(:, PMIN);
Pmax = mpc.gen(:, PMAX);
Qmin = mpc.gen(:, QMIN);
Qmax = mpc.gen(:, QMAX);
rateA = mpc.branch(:, RATE_A);

%%-----  set power flow constraints  -----
%% get bus index lists of each type of bus
[ref, pv, pq] = bustypes(bus, gen);
%% generator info
on = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?
gbus = gen(on, GEN_BUS);                %% what buses are they at?
% nb = size(bus, 1);
% ngon = size(on, 1);
% Cg = sparse(gbus, (1:ngon)', 1, nb, ngon);  %% connection matrix
%                                             %% element i, j is 1 if
%                                             %% gen on(j) at bus i is ON

%% voltages at generator buses
bus(gbus, VMAX) = gen(on, VG);
bus(gbus, VMIN) = gen(on, VG);
%% bus 33 voltage to be fixed at 0.94 pu.
bus(33, VMIN) = 0.94;              % fixing the voltage of bus 33 to 0.94
bus(33, VMAX) = 0.94;
%% generator limits (except slack)
gen(on, PG) = 0;

gen(on, QMAX) = Inf;
gen(on, QMIN) = -Inf;

%% find slack gens, remove P lims
gref = find(bus(gen(:, GEN_BUS), BUS_TYPE) == REF);
gen(gref, PMAX) = Inf;
gen(gref, PMIN) = -Inf;

%% remove line lims
mpci.branch(:, RATE_A) = 0;
[mpci.bus, mpci.gen] = deal(bus, gen);
mpc = int2ext(mpci);
result = runopf(mpc, mpopt, fname, solvedcase);
Q10=result.bus(10,3)         % net reactive power injection at bus 10
Qinjgen=result.gen(2,QG)
V33=result.bus(33, VM)                   % voltage at bus 33
%% restore original limits
varargout{1}.bus(:, VMIN) = Vmin;
varargout{1}.bus(:, VMAX) = Vmax;
varargout{1}.gen(:, PMIN) = Pmin;
varargout{1}.gen(:, PMAX) = Pmax;
varargout{1}.gen(:, QMIN) = Qmin;
varargout{1}.gen(:, QMAX) = Qmax;
varargout{1}.branch(:, RATE_A) = rateA;
