function fm_pv(flag)
% FM_PV define PV generators
%
% FM_PV(FLAG)
%     FLAG  1 -> algebraic equations
%           2 -> algebraic Jacobian matrices
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    09-Jul-2003
%Update:    25-Apr-2004
%Version:   1.0.2
%
%Modifyed 2006-01-22 by Lars.Lindgren@ts.mah.se
%
%E-mail:    fmilano@thunderbox.uwaterloo.ca
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2003 Federico Milano
% 
% This toolbox is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2.0 of the License, or
% (at your option) any later version.
% 
% This toolbox is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANDABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this toolbox; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307,
% USA.

global PV Bus DAE Settings
if ~PV.n & ~PV.pq.n, return, end

switch flag
 case 1
  
  K = 1+DAE.kg*PV.con(:,10);
  DAE.gp(PV.bus) = DAE.gp(PV.bus) - K.*PV.con(:,4);
  Qg = DAE.gq(PV.bus);
  DAE.gq(PV.bus) = 0;
  
  % check reactive power limits
  idx = [];
  
  if Settings.pv2pq 
     % Limit check improved by Lars L. 2006-01.
     %moderate performence (1 extra iteration per limit) but more correct
     %the margin to the limit must be twice the maximum error in the previus iteration
     %And to avoid bad interaction only one at a time.
     %not fully tested and could be improved in many ways
     %maybe change many buses with weak interaktion at the same time.
     %or make the conversion from PV to PQ reversible.
     if isempty(DAE.g)
        prev_err=inf;
     else
        prev_err=max(abs(DAE.g));
     end
     %idx = find(Qg+prev_err*10 < PV.con(:,7)); %old check 
    %idx = find(Qg+Settings.lftol < PV.con(:,7)); %original check
     
     [tmp,idx]=max(PV.con(:,7)-Qg(:)-prev_err*2);  
     if tmp < 0
        idx=[];
     end
  end
  
  if ~isempty(idx),
    fm_disp(strcat('Switch PV bus #', num2str(PV.con(idx,1)), ...
                   ' to PQ bus: Min Qg reached'))
    
    DAE.gp(PV.bus(idx)) = DAE.gp(PV.bus(idx)) + K(idx).*PV.con(idx,4);
    DAE.gq(PV.bus(idx)) = Qg(idx);
    
    PV.pq.con = [PV.pq.con; [PV.con(idx,:),PV.con(idx,7)]];
    PV.pq.n = PV.pq.n + length(idx);
    PV.pq.bus = [PV.pq.bus;PV.bus(idx)];
    
    idx1 = idx;
    
  else
    idx1 = [];
  end
  
  if Settings.pv2pq 
     % Limit check improved by Lars L. 2006-01.
     %moderate performence (1 extra iteration per limit) but more correct
     %the margin to the limit must be twice the maximum error in the previus iteration
     %And to avoid bad interaction only one at a time.
     %not fully tested and could be improved in many ways
     %maybe change many buses with weak interaktion at the same time.
     %or make the conversion from PV to PQ reversible.
     %idx = find(Qg-prev_err/10 > PV.con(:,6)); %old limit check
     %    idx = find(Qg-Settings.lftol > PV.con(:,6)); %original check
     [tmp,idx]=min(PV.con(:,6)-Qg(:)+prev_err*2);
     if tmp > 0
        idx=[];
     end
  end
  
  if ~isempty(idx), 
    fm_disp(strcat('Switch PV bus #', num2str(PV.con(idx,1)), ...
                   ' to PQ bus: Max Qg reached'))
    
    DAE.gp(PV.bus(idx)) = DAE.gp(PV.bus(idx)) + K(idx).*PV.con(idx,4);
    DAE.gq(PV.bus(idx)) = Qg(idx);
    
    PV.pq.con = [PV.pq.con; [PV.con(idx,:),PV.con(idx,6)]];
    PV.pq.n = PV.pq.n + length(idx);
    PV.pq.bus = [PV.pq.bus;PV.bus(idx)];
    
    idx2 = idx;
    
  else
    idx2 = [];
  end
  idx = [idx1; idx2];
  if ~isempty(idx), 
    PV.n = PV.n - length(idx); 
    PV.bus(idx) = []; 
    PV.con(idx,:) = []; 
  end
  
  if PV.pq.n
    K = 1+DAE.kg*PV.pq.con(:,10); 
    DAE.gp(PV.pq.bus) = DAE.gp(PV.pq.bus) - K.*PV.pq.con(:,4);
    DAE.gq(PV.pq.bus) = DAE.gq(PV.pq.bus) - PV.pq.con(:,end);    
  end
  
 case 2
  
  DAE.J21(PV.bus,:) = 0;
  DAE.J22(PV.bus,:) = 0;
  DAE.J12(:,PV.bus) = 0;
  DAE.J22(:,PV.bus) = 0;
  if Settings.octave
    DAE.J22(PV.bus,PV.bus) = DAE.J22(PV.bus,PV.bus) + eye(PV.n);    
  else
    DAE.J22(PV.bus,PV.bus) = DAE.J22(PV.bus,PV.bus) + speye(PV.n);
  end
  
end