function [bus, gen, branch] = opfsoln(baseMVA, bus0, gen0, branch0, ...
                                Ybus, Yf, Yt, V, Sg, lambda, ref, pv, pq, mpopt)
%OPFSOLN   Updates bus, gen, branch data structures to match opf soln.
%   [bus, gen, branch] = opfsoln(baseMVA, bus0, gen0, branch0, ...
%                              Ybus, Yf, Yt, V, Sg, lambda, ref, pv, pq, mpopt)

%   MATPOWER
%   $Id: opfsoln.m,v 1.4 2004/08/23 20:56:51 ray Exp $
%   by Ray Zimmerman, PSERC Cornell
%   Copyright (c) 1996-2004 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ for more info.

%% constants
j = sqrt(-1);
nl = size(branch0, 1);      %% number of lines

%% define named indices into data matrices
[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;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ...
    GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, ...
    RATE_C, TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST] = idx_brch;

%% initialize return values
bus     = bus0;
gen     = gen0;
branch  = branch0;

%%----- update bus voltages -----
bus(:, VM) = abs(V);
bus(:, VA) = angle(V) * 180 / pi;

%%----- update Pg and Qg for all gens -----
%% generator info
on = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?
gbus = gen(on, GEN_BUS);                %% what buses are they at?

%% copy back Pg & Vg
gen(:, PG) = zeros(size(gen, 1), 1);
gen(on, PG) = real(Sg) * baseMVA;
gen(on, VG) = bus(gbus, VM);

%% compute Qg if not passed in with Sg
if ~any(imag(Sg))
    %% This is slow in Matlab 5 ...
    % Sg = V(gbus) .* conj(Ybus(gbus, :) * V);
    %% ... so we do this instead ...
    temp = Ybus.';
    Sg = V(gbus) .* conj(temp(:, gbus).' * V);
    Qg = imag(Sg) + bus(gbus, QD) / baseMVA;    %% inj Q + local Qd
else
    Qg = imag(Sg);
end

%% update Qg for all generators
gen(:, QG) = zeros(size(gen, 1), 1);        %% zero out all Qg
gen(on, QG) = Qg * baseMVA;                 %% except for on-line generators

%%----- update/compute branch power flows -----
out = find(branch(:, BR_STATUS) == 0);      %% out-of-service branches
br = find(branch(:, BR_STATUS));            %% in-service branches
Sf = V(branch(br, F_BUS)) .* conj(Yf(br, :) * V) * baseMVA; %% complex power at "from" bus
St = V(branch(br, T_BUS)) .* conj(Yt(br, :) * V) * baseMVA; %% complex power injected at "to" bus
branch(br, [PF, QF, PT, QT]) = [real(Sf) imag(Sf) real(St) imag(St)];
branch(out, [PF, QF, PT, QT]) = zeros(length(out), 4);

%%----- update lambda's and mu's -----
%% sizes of things
nb = size(bus, 1);
nl = size(branch, 1);
npv = length(pv);
npq = length(pq);
ng = length(on);                        %% number of generators that are turned on
nbr = length(br);                       %% number of branches in service

%% initialize with all zeros
bus(:, [LAM_P, LAM_Q, MU_VMIN, MU_VMAX]) = zeros(nb, 4);
gen(:, [MU_PMIN, MU_PMAX, MU_QMIN, MU_QMAX]) = zeros(size(gen, 1), 4);
branch(:, [MU_SF, MU_ST]) = zeros(nl, 2);

%% set up indexing for lambda
i1 = 1;         i2 = nb;            %% i1:i2 - P mismatch, all buses
i3 = i2 + 1;    i4 = i2 + nb;       %% i3:i4 - Q mismatch, all buses
i5 = i4 + 1;    i6 = i4 + nb;       %% i5:i6 - Vmin, all buses
i7 = i6 + 1;    i8 = i6 + nb;       %% i7:i8 - Vmax, all buses
i9 = i8 + 1;    i10 = i8 + ng;      %% i9:i10 - Pmin, gen buses
i11 = i10 + 1;  i12 = i10 + ng;     %% i11:i12 - Pmax, gen buses
i13 = i12 + 1;  i14 = i12 + ng;     %% i13:i14 - Qmin, gen buses
i15 = i14 + 1;  i16 = i14 + ng;     %% i15:i16 - Qmax, gen buses
i17 = i16 + 1;  i18 = i16 + nbr;    %% i17:i18 - |Sf| line limit
i19 = i18 + 1;  i20 = i18 + nbr;    %% i19:i20 - |St| line limit

%% copy multipiers to bus, gen, branch
bus(:, LAM_P)       = lambda(i1:i2) / baseMVA;
bus(:, LAM_Q)       = lambda(i3:i4) / baseMVA;
bus(:, MU_VMIN)     = lambda(i5:i6);
bus(:, MU_VMAX)     = lambda(i7:i8);

gen(on, MU_PMIN)    = lambda(i9:i10) / baseMVA;
gen(on, MU_PMAX)    = lambda(i11:i12) / baseMVA;
gen(on, MU_QMIN)    = lambda(i13:i14) / baseMVA;
gen(on, MU_QMAX)    = lambda(i15:i16) / baseMVA;

branch(br, MU_SF)   = lambda(i17:i18) / baseMVA;
branch(br, MU_ST)   = lambda(i19:i20) / baseMVA;

return;
