function MargRespObj = margplot(RespObj,Systems,X,Y,PlotStyle,MarginVals);
%MARGPLOT Create a Stability Margin Response Object for the Control System Toolbox

%   Authors: Karen Gondoly, 2-13-98
%   Revised: Adam W. DiVergilio, 11-99
%   Copyright 1986-2000 The MathWorks, Inc. 
%   $Revision: 1.22 $ $Date: 2000/06/19 21:02:56 $

ni = nargin;

%---Properties
MargRespObj = struct(...
   'FrequencyUnits','',...
   'MagnitudeUnits','',...
   'PhaseUnits','',...
   'StabilityMarginValue',[]);

%---Return cases (ni==0,1)
switch ni
case 0
   RespObj = response;
   MargRespObj = class(MargRespObj,'margplot',RespObj);
   return
case 1
   if ~isa(RespObj,'response')
      error('Invalid Response Object passed to MARGPLOT')
   elseif isa(RespObj,'margplot')
      MargRespObj = RespObj;
   else
      MargRespObj = class(MargRespObj,'margplot',RespObj);
   end
   return
end

%---Initialize MargRespObj
MargFlag = isa(RespObj,'margplot');
if MargFlag
   MargRespObj = RespObj;
   RespObj = MargRespObj.response;
else
   Preferences = get(RespObj,'Preferences');
   MargRespObj.FrequencyUnits = Preferences.FrequencyUnits;
   MargRespObj.MagnitudeUnits = Preferences.MagnitudeUnits;
   MargRespObj.PhaseUnits = Preferences.PhaseUnits;
   PA = get(RespObj,'PlotAxes');
   set(PA(1:2:end-1),'XScale',Preferences.FrequencyScale,'YScale',Preferences.MagnitudeScale);
   set(PA(2:2:end),  'XScale',Preferences.FrequencyScale,'YScale','linear');
   set(get(PA(1),'Title'),'Color','k');
end

%---Access needed parent property values
PV = get(RespObj,{'PlotAxes','UIContextMenu','SystemNames'});
[LTIdisplayAxes,ContextMenu,SystemNames] = deal(PV{:});
if MargFlag
   % Adding to existing MARGPLOT object
   RH0 = get(RespObj,'ResponseHandles');
else
   % New MARGPLOT object
   RH0 = {};
end

%---Get dimension info and names of new systems
[numout,numin] = size(LTIdisplayAxes);
NumOldSys = length(RH0);                      % number of existing systems
NewSystemNames = SystemNames(NumOldSys+1:end);
NumNewSys = length(NewSystemNames);           % number of added systems

%---Store system names in MARGINVALS structure
MarginVals.System = NewSystemNames;

%---Initialize ResponseHandles array
ResponseHandles = cell(NumNewSys,1);

%---Get the style and legend
Style = struct2cell(PlotStyle);
[CC,LL,MM,LegendStr] = deal(Style{:});

%---Add legend to system name for new systems
set(ContextMenu.Systems.Names(NumOldSys+1),...
   'Label',sprintf('%s (%s)',NewSystemNames{1},LegendStr),...
   'ForegroundColor',CC)

%---Convert to preferred units
X = unitconv(X{1}{1},'rad/sec',MargRespObj.FrequencyUnits);
Magdata = unitconv(Y{1},'abs',MargRespObj.MagnitudeUnits);
Phdata = unitconv(Y{2},'deg',MargRespObj.PhaseUnits);
MarginVals.GMFrequency = unitconv(MarginVals.GMFrequency,'rad/sec',MargRespObj.FrequencyUnits);
MarginVals.PMFrequency = unitconv(MarginVals.PMFrequency,'rad/sec',MargRespObj.FrequencyUnits);
zeroGM = (MarginVals.GainMargin == 0);
MarginVals.GainMargin  = unitconv(MarginVals.GainMargin, 'db',MargRespObj.MagnitudeUnits);
zeroPM = (MarginVals.PhaseMargin == 0);
MarginVals.PhaseMargin = unitconv(MarginVals.PhaseMargin,'deg',MargRespObj.PhaseUnits);
p180 = unitconv(180,'deg',MargRespObj.PhaseUnits);
m0 = unitconv(0,'db',MargRespObj.MagnitudeUnits);

switch lower(MargRespObj.MagnitudeUnits)
case {'db','decibels'}
   GmPlot = -1*MarginVals.GainMargin;
otherwise
   GmPlot = 1/MarginVals.GainMargin;
end

%---Determine if system is discrete
dflag = abs(Systems{1}.Ts);

%---Initializations
RH1 = cell(1);
RH2 = cell(1);
MagLine = 1;
PhLine = 1;

%---Response line UserData
LineUd = struct('System',NewSystemNames{1},...
   'Input',1,'Output',1,'Array',[1 1]);

%---Magnitude
if ndims(Magdata)>2
   Magdata = squeeze(Magdata(1,1,:));
end

%---Magnitude
RH1{1}(1,1)=line(X,Magdata,...
   'Color',CC,'Marker',MM,'LineStyle',LL,'UserData',LineUd,...
   'Tag','LTIresponseLines','DeleteFcn','delresp(gcbo)',...
   'Visible','on','Parent',LTIdisplayAxes(1,1),...
   'ButtonDownFcn','rguifcn(''showbox'',gcbf);');

%---Phase
if ndims(Phdata)>2
   Phdata = squeeze(Phdata(1,1,:));
end
RH2{1}(1,1)=line(X,Phdata,...
   'Color',CC,'Marker',MM,'LineStyle',LL,'UserData',LineUd,...
   'Tag','LTIresponseLines','DeleteFcn','delresp(gcbo)',...
   'Visible','on','Parent',LTIdisplayAxes(2,1),...
   'ButtonDownFcn','rguifcn(''showbox'',gcbf);');

%---Nyquist line
if dflag>0
   %---This is currently never called since the system is not passed to margin
   % Only draw Nyquist line when data extends, at least, up to NyqFreq
   NyqFreq = unitconv(pi/dflag,'rad/sec',MargRespObj.FrequencyUnits);
   if abs(X(end)-NyqFreq)<NyqFreq/10
      RH1{1}(MagLine+1,1)=line([NyqFreq;NyqFreq],[0;1],'Color','k','Linestyle','-',...
         'Parent',LTIdisplayAxes(1,1),'HandleVis','off','HitTest','off','Tag','NyquistLines');
      MagLine = MagLine+1;
      RH2{1}(PhLine+1,1)=line([NyqFreq;NyqFreq],[0;1],'Color','k','Linestyle','-',...
         'Parent',LTIdisplayAxes(2,1),'HandleVis','off','HitTest','off','Tag','NyquistLines');
      PhLine = PhLine+1;
   end
end

% Determine -180 (mod 360) phase lines associated gain and phase margins
Xlim = [min(X) max(X)];
Wcg = MarginVals.GMFrequency;
Wcp = MarginVals.PMFrequency;
if Wcg>=Xlim(1) & Wcg<=Xlim(2)
   GmPhaseLine = p180 * round(interp1(X(:),Phdata(:),Wcg)/p180);
end
if Wcp>=Xlim(1) & Wcp<=Xlim(2)
   PmPhaseLine = p180 * round((interp1(X(:),Phdata(:),Wcp)-MarginVals.PhaseMargin)/p180);
end

%---Gm lines
if Xlim(1)<=Wcg & Xlim(2)>=Wcg & ~zeroGM
   RH1{1}(MagLine+1,1)=line([Wcg;Wcg],[GmPlot;m0],...
      'Parent',LTIdisplayAxes(1),'Tag','GainMarginLine',...
      'HandleVisibility','off','HitTest','off',...
      'Color','k','LineStyle','-');
   RH1{1}(MagLine+2,1)=line(Xlim,[m0;m0],...
      'Parent',LTIdisplayAxes(1),'Tag','GainMarginHorizontalLine',...
      'HandleVisibility','off','HitTest','off',...
      'Color','k','LineStyle',':');
   RH1{1}(MagLine+3,1)=line([Wcg;Wcg],[m0;m0],...
      'Parent',LTIdisplayAxes(1),'Tag','GainMarginVerticalLine',...
      'HandleVisibility','off','HitTest','off',...
      'Color','k','LineStyle',':');
   RH2{1}(PhLine+1,1)=line([Wcg;Wcg],[GmPhaseLine;GmPhaseLine],...
      'Parent',LTIdisplayAxes(2),'Tag','GainMarginVerticalLine',...
      'HandleVisibility','off','HitTest','off',...
      'Color','k','LineStyle',':');
   RH2{1}(PhLine+2,1)=line(Xlim,[GmPhaseLine;GmPhaseLine],...
      'Parent',LTIdisplayAxes(2),'Tag','GainMarginHorizontalLine',...
      'HandleVisibility','off','HitTest','off',...
      'Color','k','LineStyle',':');
   MagLine = MagLine+3;
   PhLine = PhLine+2;
end

%---Pm lines
if Xlim(1)<=Wcp & Xlim(2)>=Wcp & ~zeroPM
   RH2{1}(PhLine+1,1)=line([Wcp;Wcp],[PmPhaseLine;MarginVals.PhaseMargin+PmPhaseLine],...
      'Parent',LTIdisplayAxes(2),'Tag','PhaseMarginLine',...
      'HandleVisibility','off','HitTest','off',...
      'Color','k','LineStyle','-');
   RH2{1}(PhLine+2,1)=line(Xlim,[PmPhaseLine;PmPhaseLine],...
      'Parent',LTIdisplayAxes(2),'Tag','PhaseMarginHorizontalLine',...
      'HandleVisibility','off','HitTest','off',...
      'Color','k','LineStyle',':');
   RH2{1}(PhLine+3,1)=line([Wcp;Wcp],[PmPhaseLine;PmPhaseLine],...
      'Parent',LTIdisplayAxes(2),'Tag','PhaseMarginVerticalLine',...
      'HandleVisibility','off','HitTest','off',...
      'Color','k','LineStyle',':');
   RH1{1}(MagLine+1,1)=line([Wcp;Wcp],[m0;m0],...
      'Parent',LTIdisplayAxes(1),'Tag','PhaseMarginVerticalLine',...
      'HandleVisibility','off','HitTest','off',...
      'Color','k','LineStyle',':');
   RH1{1}(MagLine+2,1)=line(Xlim,[m0;m0],...
      'Parent',LTIdisplayAxes(1),'Tag','PhaseMarginHorizontalLine',...
      'HandleVisibility','off','HitTest','off',...
      'Color','k','LineStyle',':');
end

ResponseHandles{1} = {RH1;RH2};

%---Gain text
if strncmpi(MargRespObj.MagnitudeUnits,'d',1)
   mstr = sprintf(' %s',MargRespObj.MagnitudeUnits);
else
   mstr = sprintf('');
end
if isinf(MarginVals.GainMargin)
   tStr = sprintf('Gm = Inf');
elseif ~unitconv(MarginVals.GainMargin,MargRespObj.MagnitudeUnits,'db')
   tStr = sprintf('Gm = %0.5g%s',MarginVals.GainMargin,mstr);
else
   tStr = sprintf('Gm = %0.5g%s (at %0.5g %s)',...
      MarginVals.GainMargin,mstr,...
      MarginVals.GMFrequency,MargRespObj.FrequencyUnits);
end

%---Phase text
if isinf(MarginVals.PhaseMargin)
   tStr = sprintf('%s,  Pm = Inf',tStr);
elseif ~MarginVals.PhaseMargin
   tStr = sprintf('%s,  Pm = 0 (unstable closed loop)',tStr);
else
   tStr = sprintf('%s,  Pm = %0.5g %s (at %0.5g %s)',tStr,...
      MarginVals.PhaseMargin,MargRespObj.PhaseUnits,...
      Wcp,MargRespObj.FrequencyUnits);
end

%---Display title string
set(get(LTIdisplayAxes(1,1),'Title'),'String',tStr);

%---Set the BackgroundAxes labels
BackAxLabels = {...
   'Title','Bode Diagram',...
   'Ylabel',sprintf('Phase (%s); Magnitude (%s)',MargRespObj.PhaseUnits,MargRespObj.MagnitudeUnits),...
   'Xlabel',sprintf('Frequency (%s)',MargRespObj.FrequencyUnits)};

if MargFlag  % Append Data
   RespObj = pvset(RespObj,'UIContextMenu',ContextMenu,...
      'ResponseHandles',[RH0;ResponseHandles],...
      BackAxLabels{:});
   MargRespObj.StabilityMarginValue = [MargRespObj.StabilityMarginValue;MarginVals];
   MargRespObj.response = RespObj;

else
   %---Hide irrelevant menus
   set([ContextMenu.PlotOptions.Main,ContextMenu.ChannelMenu,...
         ContextMenu.GroupMenu.Main],'visible','off');
   set(ContextMenu.GridMenu,'Separator','on')
   RespObj = pvset(RespObj,'ResponseHandles',ResponseHandles,...
      'UIContextMenu',ContextMenu,BackAxLabels{:});
   MargRespObj.StabilityMarginValue = MarginVals;
   MargRespObj = class(MargRespObj,'margplot',RespObj);
end

%---Update HG pointer to Response Object
builtin('set',ContextMenu.Main,'userdata',MargRespObj)
