function drawPVcurves(casename, loadvarloc, corrected_list, combined_list, flag_combinedCurve, busesToDraw)
% function: draw PV curves for specified buses
% [INPUT PARAMETERS]
% corrected_list, combined_list: data points obtained from CPF solver
% loadvarloc: load variation location(in external bus numbering). Single bus supported so far.
% flag_combinedCurve: flag indicating if the prediction-correction curve will be drawn
% busesToDraw: bus indices whose PV curve will be be drawn
% created by Rui Bo on 2008/01/13

%% assign default parameters
if nargin < 6
    busesToDraw = loadvarloc;   % draw the curve for the load changing bus
end
if isempty(busesToDraw)
    busesToDraw = loadvarloc;   % draw the curve for the load changing bus
end

%% load the case & convert to internal bus numbering
[baseMVA, bus, gen, branch] = feval(casename);
nb = size(bus, 1);

correctedDataNum = size(corrected_list, 2) - 1;
combinedDataNum  = size(combined_list, 2) - 1;

%% prepare data for drawing
lambda_corrected = corrected_list(nb+1, [2:correctedDataNum+1]);
lambda_combined  = combined_list(nb+1, [2:combinedDataNum+1]);

fprintf('Start plotting CPF curve(s)...\n');
for j = 1:length(busesToDraw)%for i = 1+npv+1:1+npv+npq
    i = find(corrected_list(:, 1) == busesToDraw(j)); % find row index
    
    %% get voltage magnitudes
    Vm_corrected = abs(corrected_list(i, [2:correctedDataNum+1]));
    Vm_combined  = abs(combined_list(i, [2:combinedDataNum+1]));
    
    %% create a new figure
    figure;
    hold on;
    
    %% plot PV curve
    plot(lambda_corrected, Vm_corrected, 'bx-');
    
    %% plot CPF prediction-correction curve
    if flag_combinedCurve == true
        plot(lambda_combined, Vm_combined, 'r.-');
        legend('CPF Curve', 'Prediction-Correction Curve');
        legend('Location', 'Best');
    end
    
    %% add plot title
    title(['Vm at bus ' int2str(busesToDraw(j)) ' w.r.t. load (p.u.) at ' int2str(loadvarloc)]);
end
fprintf('Plotting is done.\n');
