function dummy = refreshpz(Editor,event,varargin)
%REFRESHPZ  Refreshes plot after compensator update.

%   Author(s): P. Gahinet
%   Copyright 1986-2000 The MathWorks, Inc. 
%   $Revision: 1.14 $  $Date: 2000/08/17 22:16:18 $

% REVISIT
dummy = [];

persistent PZListener 

% Process events
switch event 
case 'init'
   % Initialization for dynamic gain update (drag).
   PZGroup = varargin{1};     % Modified PZGROUP
   
   % Find related PZVIEW objects
   HG = Editor.HG;
   ig = find(Editor.LoopData.Compensator.PZGroup==PZGroup); 
   PZView(1) = HG.Compensator.Magnitude(ig);
   PZView(2) = HG.Compensator.Phase(ig);
   
   % Switch editor's RefreshMode to quick
   Editor.RefreshMode = 'quick';
   
   % Save initial data (units = rad/sec,abs,deg)
   InitData = struct(...
       'PZGroup',get(PZGroup),...
       'Frequency',Editor.Frequency,...
       'Magnitude',Editor.Magnitude,...
       'Phase',Editor.Phase);
   
   % Install listener on PZGROUP data
   PZListener = handle.listener(PZGroup,'PZDataChanged',...
       {@LocalUpdatePlot Editor InitData PZView});
   
case 'finish'
   % Clean up after dynamic gain update (drag)
   % Return editor's RefreshMode to normal
   Editor.RefreshMode = 'normal';
   
   % Delete listener
   delete(PZListener);
   PZListener = [];
end


%-------------------------Local Functions-------------------------

%%%%%%%%%%%%%%%%%%%%%%%
%%% LocalUpdatePlot %%%
%%%%%%%%%%%%%%%%%%%%%%%
function LocalUpdatePlot(PZGroup,event,Editor,InitData,PZView)
% Update plot

% RE:  * PZGroup: current PZGROUP data
%      * Working units are (rad/sec,abs,deg)

HG = Editor.HG;
LoopData = Editor.LoopData;
Ts = getTs(LoopData);
Format = getformat(LoopData);

% Conversion factor
FreqConvert = unitconv(1,'rad/sec',Editor.FrequencyUnits);

% Construct corrective factor
zcorr = [InitData.PZGroup.Pole ; PZGroup.Zero];
pcorr = [InitData.PZGroup.Zero ; PZGroup.Pole];

% Natural and peaking frequencies for new pole/zero locations (in rad/sec)
[W0,Zeta] = damp([PZGroup.Zero ; PZGroup.Pole],Ts);
if Ts,
	% Keep root freq. below Nyquist freq.
	W0 = min(W0,pi/Ts);
end
t = W0.^2 .* (1 - 2 * Zeta.^2);
Wpeak = sqrt(t(t>0,:));

% Update mag and phase data
% RE: Update editor properties (used by INTERPY and REFRESHMARGIN)
GainMag = abs(getgain(LoopData));
InitMag = [InitData.Magnitude ; ...
        interpmag(Editor,InitData.Frequency,InitData.Magnitude,[Wpeak;W0])];
InitPhase = [InitData.Phase ; interp1(InitData.Frequency,InitData.Phase,[Wpeak;W0])];
[W,iu] = LocalUniqueWithinTol([InitData.Frequency;Wpeak;W0],1e3*eps);  % sort + unique
[Mag,Phase] = LocalBodeUpdate(W,InitMag(iu),InitPhase(iu),zcorr,pcorr,Ts,Format);
Editor.Frequency = W;
Editor.Magnitude = Mag;
Editor.Phase = Phase;

% Update plot
MagData = unitconv(GainMag*Editor.Magnitude,'abs',Editor.MagnitudeUnits);
set(HG.BodePlot(1),'Xdata',FreqConvert*Editor.Frequency,'Ydata',MagData)
set(HG.BodePlot(2),'Xdata',FreqConvert*Editor.Frequency,...
    'Ydata',unitconv(Editor.Phase,'deg',Editor.PhaseUnits))

% Update X location of moved roots
% REVISIT: merge
W0 = FreqConvert * W0;
p = PZView(1);  % mag
set([p.Zero;p.Pole],{'Xdata'},num2cell(W0,2));
p = PZView(2);  % phase
set([p.Zero;p.Pole],{'Xdata'},num2cell(W0,2));

% Update location of notch width markers
if strcmp(PZGroup.Type,'Notch')
    Wm = FreqConvert*notchwidth(PZGroup,Ts);
    % Markers
    Extras = get(PZView(1),'Extra');
    set(Extras,{'Xdata'},num2cell(Wm,2));
    % Rulers
    Rulers = get(PZView(1),'Ruler');
    set(Rulers,{'Xdata'},{ Wm([1 1],:) ; Wm([2 2],:) });
end

% Interpolate Y values
PhaseData = unitconv(Editor.Phase,'deg',Editor.PhaseUnits);
interpy(Editor,MagData,PhaseData);

% Update stability margins (using interpolation)
refreshmargin(Editor);



%%%%%%%%%%%%%%%%%%%
% LocalBodeUpdate %
%%%%%%%%%%%%%%%%%%%
function [mag,phase] = LocalBodeUpdate(w,mag,phase,zcorr,pcorr,Ts,Format)
%BODEUPDATE  Updates mag and phase data when moving poles and zeros.
%
%   [MAG,PHASE] = BODEUPDATE(W,MAG,PHASE,ZCORR,PCORR,TS,FORMAT)
% 
%   Input:  
%     * MAG,PHASE  :  initial data
%     * ZCORR,PCORR:  corrective term in ZP format
%
%   Returns the updated mag and phase when the inital transfer 
%   function is multiplied by 
%     * prod(s-zj)/prod(s-pj)       if Format = 'ZeroPoleGain'
%     * prod(1-s/zj)/prod(1-s/pj)   if Format = 'TimeConstant'.

% S or Z vectors
s = i*w(:).';
if Ts,
   s = exp(Ts * s);
end
ls = length(s);

% Compute correction
switch lower(Format(1))
case 'z'  % zero/pole/gain
   % Corrective term is prod(s-zj)/prod(s-pj)
   sz = s(ones(1,length(zcorr)),:) - zcorr(:,ones(1,ls));
   sp = s(ones(1,length(pcorr)),:) - pcorr(:,ones(1,ls));
   a = prod(sz,1);
   b = prod(sp,1);
   Correction = ones(size(a));
   nzb = find(b);
   Correction(:,nzb) = a(:,nzb)./b(:,nzb);
   
case 't'  % time constant
   % Corrective term is 
   %    s^m * prod(1-s/zj)/prod(1-s/pj)  in continuous time
   %    w^m * prod(1-w/(zj-1))/prod(1-w/(pj-1))  in discrete time with w=z-1
   if Ts
      zcorr = zcorr-1;   pcorr = pcorr-1;   s = s-1;
   end
   lz = length(zcorr);
   zcorr(~zcorr,:) = [];
   lp = length(pcorr);
   pcorr(~pcorr,:) = [];
   sz = 1 - s(ones(1,length(zcorr)),:) ./ zcorr(:,ones(1,ls));
   sp = 1 - s(ones(1,length(pcorr)),:) ./ pcorr(:,ones(1,ls));
   m = lz - length(zcorr) - lp + length(pcorr);
   
   % Guard against divide by zero (pole/zero at s=0)
   a = prod(sz,1);
   b = prod(sp,1);
   Correction = ones(size(a));
   nzb = find(b);
   Correction(:,nzb) = s(:,nzb).^m .* a(:,nzb)./b(:,nzb);
   
end

% Update mag and phase
mag = mag .* abs(Correction(:));
phase = phase + (180/pi) * unwrap(angle(Correction(:)));


%%%%%%%%%%%%%%%%%%%%%%%%
% LocalUniqueWithinTol %
%%%%%%%%%%%%%%%%%%%%%%%%
function [w,iu] = LocalUniqueWithinTol(w,rtol)
% Eliminates duplicates within RTOL (relative tolerance)
% Helps prevent reintroducing duplicates during unit conversions

% Sort W
[w,iu] = sort(w);

% Eliminate duplicates
lw = length(w);
dupes = find(w(2:lw)-w(1:lw-1)<=rtol*w(2:lw));
w(dupes,:) = [];
iu(dupes,:) = [];
