www.pudn.com > Closid30.rar > CLPLOTEV.M, change:1998-10-22,size:28594b


function clplotev(option,eval) 
 
%Draws figures of the evaluation methods. 
%Input:option==1 closed-loop transfer functions 
%            ==2 closed-loop poles 
%            ==3 input-output simulation 
%            ==4 correlation tests 
%            ==5 step responses 
%            ==6 plant transfer  
%            ==7 plant pole-zeros 
%        eval=='all' all models are calculated and plotted 
%            =='one' the validation result of one model is 
%                      added or deleted from the evaluation figure 
% 
% 09-06-1997 
% (c) Edwin van Donkelaar, Paul Van den Hof 
% Mechanical Engineering Systems and Control Group 
% Delft University of Technology 
% Last update: 19-05-1998 
%              24-08-1998 Correction closed-loop f-response plots for  
%                         incorporation of sampling time unequal to 1. 
%              25-08-1998 Correction correlation test on u. 
%              15-09-1998 Adapted to handle validation data. 
%              21-10-1998 Removing th2ss in closed-loop pole calculation, 
%                         and replacing by polynomial calculation. 
% 
%********************************************************************** 
 
global CLIDmods CLIDscrz CLIDevas 
 
z=CLIDscrz; 
 
% Colors 
[frmcol,txtcol,edicol,pshcol,radcol,axscol,c7,c8,c9,c10,... 
modcol,edtcol,pubcol,putcol,chbcol,chtcol,txbcol,pstcol]=clsetcol; 
 
% Total number of models, active models and indices  
% find all active models 
allmodels=CLIDmods(7,:); 
modacta=find(CLIDmods(7,:)==1);   
% find active parametric models 
modact=find(CLIDmods(57,modacta)==1);  
modact=modacta(modact); 
totmod=length(modact); 
% find active frequency domain models 
modfact=find(CLIDmods(57,modacta)==0);  
modfact=modacta(modfact); 
totfmod=length(modfact); 
 
% Total number of models, active and indices  
allmodels=CLIDmods(7,:); 
totmod=sum(CLIDmods(7,:)); 
modact=find(CLIDmods(7,:)==1); 
 
% Find total number of active controllers and their indices 
allcons=CLIDmods(17,:); 
totcon=sum(allcons); 
conact=find(CLIDmods(17,:)==1);   
 
% If the number of controllers is larger than the number of models,  
% the collors of the controllers are used in the evaluation plots. 
if totcon>totmod 
  mcact=conact; 
else 
  mcact=modact;   
end 
 
if totmod==0 
  clerrdia(['ERROR: no active model available'],1) 
  return 
end 
 
% Extract necessary data from the user data of the figure 
fignum=findobj('name',CLIDevas{option}); 
 
% Read active controller for all methods that require this  
% information. 
if option==1|option==2|option==3|option==4|option==5 
  if totcon==0 
   clerrdia(['ERROR: no active controller available'],1) 
   return 
  end 
end 
 
%********************************************************************* 
% 
%                 Closed-loop transfer functions 
% 
%********************************************************************* 
 
if option==1 
  if strcmp(eval,'all') 
    fignum=findobj(0,'name',CLIDevas{option}); 
    figure(fignum) 
    us=get(fignum,'userdata');Nus=length(us);   
 
% Get figure options and set the figure layout  
 
% Get frequency axis and channel  from user data from the figure  
    wp=us(2:Nus,1);chan=us(1,1);                    
% Get amplitude scale  
    if strcmp(get(findobj(gcf,'label','Log amplitude scale'),'checked'),'on')         
      yscale='log'; 
    else 
      yscale='lin'; 
    end 
% Get amplitude scale  
%    if us(1,2)==1 
%	  yscale='log'; 
%	else 
%	  yscale='lin'; 
%	end 
% Get grid option 
    if strcmp(get(findobj(gcf,'label','Grid'),'checked'),'on')         
      xgrid='on';ygrid='on'; 
    else 
      xgrid='off';ygrid='off'; 
    end 
%    if us(3,2)==1         
%      xgrid='on';ygrid='on'; 
%    else 
%      xgrid='off';ygrid='off'; 
%    end 
 
% Clear the axes and set the figure layout  
    if chan==1|chan==2|chan==3|chan==4|chan==7|chan==8 
      subplot(211),cla,set(gca,'yscale',yscale,'xscale','log','xgrid',xgrid,'ygrid',ygrid),ylabel('amplitude'),xlabel('frequency (rad/sec)'),hold on 
      subplot(212),cla,set(gca,'yscale','lin','xscale','log','xgrid',xgrid,'ygrid',ygrid),ylabel('phase'),xlabel('frequency (rad/sec)'),hold on          
    elseif chan==5 
      subplot(221),cla,set(gca,'yscale',yscale,'xscale','log','xgrid',xgrid,'ygrid',ygrid),ylabel('amplitude'),hold on 
      subplot(222),cla,set(gca,'yscale',yscale,'xscale','log','xgrid',xgrid,'ygrid',ygrid),hold on 
      subplot(223),cla,set(gca,'yscale',yscale,'xscale','log','xgrid',xgrid,'ygrid',ygrid),hold on 
      xlabel('frequency (rad/sec)'),ylabel('amplitude') 
      subplot(224),cla,set(gca,'yscale',yscale,'xscale','log','xgrid',xgrid,'ygrid',ygrid),xlabel('frequency (rad/sec)'),hold on    
    elseif chan==6 
      subplot(221),cla,set(gca,'yscale','lin','xscale','log','xgrid',xgrid,'ygrid',ygrid),ylabel('phase'),hold on 
      subplot(222),cla,set(gca,'yscale','lin','xscale','log','xgrid',xgrid,'ygrid',ygrid),hold on 
      subplot(223),cla,set(gca,'yscale','lin','xscale','log','xgrid',xgrid,'ygrid',ygrid),hold on,xlabel('frequency (rad/sec)'),ylabel('phase') 
      subplot(224),cla,set(gca,'yscale','lin','xscale','log','xgrid',xgrid,'ygrid',ygrid),xlabel('frequency (rad/sec)'),hold on    
    end 
 
% Calculate closed-loop transfer functions for all selected models and controllers 
    [dat,tit,info,p]=clobjget('datv');tsamp=p(2,1); % Correction 24-8-98: account for sampling times other than 0. 
	for j=1:totcon 
      C=get(CLIDmods(21,conact(j)),'userdata'); 
      for i=1:totmod 
% Parametric models  
        if CLIDmods(57,modact(i))==1 
          w=wp; 
          th=get(CLIDmods(11,modact(i)),'userdata'); 
% Corrected 22-10-1998: 
%         [a,b,c,d]=th2ss(th); P=[a b;c d];  
          [numg,deng]=th2tf(th,1);[ag,bg,cg,dg]=tf2ss(numg,deng); 
          [ag,bg,cg,dg]=minreal(ag,bg,cg,dg); 
          if isempty(ag), ag=0;bg=0;cg=0; end 
          P=[ag bg;cg dg]; 
%%%%%%%%%%%%%%%%%%%%%%%%%		   
          [TPC,nTPC]=cltpc(P,C,length(P)-1,length(C)-1); 
          [a,b,c,d]=clsplit(TPC,nTPC); 
          [a,b,c,d]=minreal(a,b,c,d); 
          [m1,f1]=dbode(a,b,c,d,tsamp,1,w);  % Correction 24-8-98: 1 -> tsamp 
          [m2,f2]=dbode(a,b,c,d,tsamp,2,w);  % Correction 24-8-98: 1 -> tsamp 
           if chan==7|chan==8 
% Corrected 22-10-1998: 
%            [ah,bh,ch,dh,kh]=th2ss(th); 
%            [ah,bh,ch,dh]=minreal(ah,kh,ch,1); 
             [numh,denh]=th2tf(th,-1);[ah,bh,ch,dh]=tf2ss(numh,denh); 
             [ah,bh,ch,dh]=minreal(ah,bh,ch,dh); 
             if isempty(ah), ah=0;bh=0;ch=0; end 
%%%%%%%%%%%%%%%%%%%%%%%%%%              
% e1 -> u: -HCS 
            [ah1,bh1,ch1,dh1]=series(ah,bh,ch,dh,a,b(:,1),-c(2,:),-d(2,1)); 
            [mh1,fh1]=dbode(ah1,bh1,ch1,dh1,tsamp,1,w);    % Correction 24-8-98: 1 -> tsamp 
% e1 -> y: HS 
            [ah2,bh2,ch2,dh2]=series(ah,bh,ch,dh,a,b(:,2),c(2,:),d(2,2)); 
            [mh2,fh2]=dbode(ah2,bh2,ch2,dh2,tsamp,1,w);    % Correction 24-8-98: 1 -> tsamp 
          end  
% nonparametric models 
        elseif CLIDmods(57,modact(i))==0 
          th=get(CLIDmods(11,modact(i)),'userdata'); 
          [w,m,f]=getff(th); 
          g=m.*exp(sqrt(-1)*f*pi/180);      
          [ac,bc,cc,dc]=clsplit(C);[mc,fc]=dbode(ac,bc,cc,dc,1,tsamp,w);  % Correction 24-8-98: 1 -> tsamp 
          gc=mc.*exp(sqrt(-1)*fc*pi/180); 
          gs=1./(ones(size(w))+gc.*g ); 
          m1=[abs(g.*gc.*gs) abs(gc.*gs)];f1=[180/pi*unwrap(angle(g.*gc.*gs))  180/pi*unwrap(angle(gc.*gs))];      
          m2=[abs(g.*gs) abs(gs)];f2=[180/pi*unwrap(angle(g.*gs)) 180/pi*unwrap(angle(gc.*gs))];      
          mh1=m1(:,2);fh1=f1(:,2)-180; 
          mh2=m2(:,2);fh2=f2(:,2); 
        end 
% Plot the models according to the selected channel 
% r1 -> u: S  
        if chan==1 
          subplot(211),plot(w,m2(:,2)),title('S: r1 -> u') 
%          if yscale=='log', axis([min(w)/5 max(w)*5 1e-3 10]),  
%		  end 
%          axis([min(w)/10 max(w)*10 1e-4 1e4]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
          subplot(212),plot(w,f2(:,2)) 
%          axis([min(w)/10 max(w)*10 -1e3 1e3]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
% r1 -> y: PS  
        elseif chan==2 
          subplot(211),plot(w,m2(:,1)),title('PS: r1 -> y') 
%          axis([min(w)/10 max(w)*10 1e-4 1e4]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
          subplot(212),plot(w,f2(:,1)) 
%          axis([min(w)/10 max(w)*10 -1e3 1e3]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
% r2 -> u: CS  
        elseif chan==3 
          subplot(211),plot(w,m1(:,2)),title('CS: r2 -> u') 
%          axis([min(w)/10 max(w)*10 1e-4 1e4]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
          subplot(212),plot(w,f1(:,2)) 
%          axis([min(w)/10 max(w)*10 -1e3 1e3]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
% r2 -> y: T  
        elseif chan==4 
          subplot(211),plot(w,m1(:,1)),title('T: r2 -> y') 
%		  if yscale=='log', axis([min(w)/5 max(w)*5 1e-3 10]),  
%		  end 
%          axis([min(w)/10 max(w)*10 1e-4 1e4]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
          subplot(212),plot(w,f1(:,1)) 
%          axis([min(w)/10 max(w)*10 -1e3 1e3]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
%[r1 r2] -> [y u] amplitude 
        elseif chan==5 
          subplot(221),plot(w,m1(:,1)),title('T: r2 -> y') 
          if yscale=='log', axis([min(w)/5 max(w)*5 1e-3 10]),  
		  end 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
          subplot(222),plot(w,m2(:,1)),title('GS: r1 -> y') 
          if yscale=='log', axis([min(w)/5 max(w)*5 1e-3 10]),  
		  end 
%          axis([min(w)/10 max(w)*10 1e-3 1e2]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
          subplot(223),plot(w,m1(:,2)),title('CS: r2 -> u') 
          if yscale=='log', axis([min(w)/5 max(w)*5 1e-3 10]),  
		  end 
%          axis([min(w)/10 max(w)*10 1e-3 1e2]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
          subplot(224),plot(w,m2(:,2)),title('S: r1 -> u') 
          if yscale=='log', axis([min(w)/5 max(w)*5 1e-3 10]),  
		  end 
%          axis([min(w)/10 max(w)*10 1e-3 1e2]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
%[r1 r2] -> [y u] phase 
        elseif chan==6  
          subplot(221),semilogx(w,f1(:,1)),title('T: r2 -> y') 
          axis([min(w)/5 max(w)*5 1e-3 100]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
          subplot(222),semilogx(w,f2(:,1)),title('GS: r1 -> y') 
          axis([min(w)/5 max(w)*5 1e-3 100]) 
%          axis([min(w)/10 max(w)*10 -1e3 1e3]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
          subplot(223),semilogx(w,f1(:,2)),title('CS: r2 -> u') 
		  axis([min(w)/5 max(w)*5 1e-3 100]) 
%          axis([min(w)/10 max(w)*10 -1e3 1e3]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
          subplot(224),semilogx(w,f2(:,2)),title('S: r1 -> u') 
          axis([min(w)/5 max(w)*5 1e-3 100]) 
%          axis([min(w)/10 max(w)*10 -1e3 1e3]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
% e1 -> u: -HCS  
        elseif chan==7 
          subplot(211),loglog(w,mh1),title('-HCS: e1 -> u') 
%          axis([min(w)/10 max(w)*10 1e-4 1e4]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
          subplot(212),semilogx(w,fh1) 
%          axis([min(w)/10 max(w)*10 -1e3 1e3]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
% e1 -> y: HS  
        elseif chan==8 
          subplot(211),loglog(w,mh2),title('HS: e1 -> y') 
%          axis([min(w)/10 max(w)*10 1e-4 1e4]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on 
          subplot(212),semilogx(w,fh2) 
%          axis([min(w)/10 max(w)*10 -1e3 1e3]) 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:),'visible','on'); 
          hold on      
        end 
      end 
    end 
% Put us in userdata of the figure   
    set(fignum,'userdata',us) 
    set(fignum,'visible','on')                           
% Get zoom option 
    if strcmp(get(findobj(gcf,'label','Zoom'),'checked'),'on')         
     zoom on 
    else 
     zoom off 
    end 
% Get zoom option 
%    if us(2,2)==1         
%      zoom on 
%    else 
%      zoom off 
%    end 
  end 
%********************************************************************* 
% 
%                 Closed-loop poles 
% 
%********************************************************************* 
elseif option==2 
  if strcmp(eval,'all') 
    fignum=findobj(0,'name',CLIDevas{option}); 
    figure(fignum),cla 
    us=get(fignum,'userdata'); 
    chan=us(1,1); 
% 
%    set(fignum,'Color',[1 1 1],'Box','on') 
 
% Set figure layout right 
% Get zoom option 
    if strcmp(get(findobj(gcf,'label','Zoom'),'checked'),'on')         
      zoom on 
    else 
      zoom off 
    end 
% Get grid option 
    if strcmp(get(findobj(gcf,'label','Grid'),'checked'),'on')         
      grid on 
    else 
      grid off 
    end 
 
% Calculate and plot pole-zeros plots for all models and controllers 
    for j=1:totcon 
      C=get(CLIDmods(21,conact(j)),'userdata'); 
      for i=1:totmod 
% Only for parametric models 
        if CLIDmods(57,modact(i))==1 
          th=get(CLIDmods(11,modact(i)),'userdata'); 
% Corrected 21-10-1998: 
          [numg,deng]=th2tf(th,1);[ag,bg,cg,dg]=tf2ss(numg,deng); 
          [ag,bg,cg,dg]=minreal(ag,bg,cg,dg); 
          if isempty(ag), ag=0;bg=0;cg=0; end 
%         [a,b,c,d]=th2ss(th); P=[a b;c d]; corrected 21-10-1998 
          P=[ag bg;cg dg]; 
          [TPC,nTPC]=cltpc(P,C,length(P)-1,length(C)-1); 
          [a,b,c,d]=clsplit(TPC,nTPC); 
          [a,b,c,d]=minreal(a,b,c,d); 
% r1 -> u: S  
          if chan==1 
% Addition 22-10-1998 (as test): 
            a1=a;b1=b(:,2);c1=c(2,:);d1=d(2,2); 
			[a1,b1,c1,d1]=minreal(a1,b1,c1,d1); 
			pols=eig(a1); 
            zers=tzero(a1,b1,c1,d1);			 
%            pols=eig(a); 
%            zers=tzero(a,b(:,2),c(2,:),d(2,2)); 
% r1 -> y: PS  
          elseif chan==2 
            pols=eig(a); 
            zers=tzero(a,b(:,2),c(1,:),d(1,2)); 
% r2 -> u: CS  
          elseif chan==3 
            pols=eig(a); 
            zers=tzero(a,b(:,1),c(2,:),d(2,1)); 
% r2 -> y: CPS  
          elseif chan==4 
            pols=eig(a); 
            zers=tzero(a,b(:,1),c(1,:),d(1,1)); 
% [r2 r1] -> [y u] 
          elseif chan==5 
            pols=eig(a); 
            zers=tzero(a,b,c,d); 
% e1 -> u: -HCS  
          elseif chan==7 
%            [ah,bh,ch,dh,kh]=th2ss(th);         
%            [ah,bh,ch,dh]=minreal(ah,kh,ch,1); corrected 21-10-1998: 
             [numh,denh]=th2tf(th,-1);[ah,bh,ch,dh]=tf2ss(numh,denh); 
             [ah,bh,ch,dh]=minreal(ah,bh,ch,dh); 
             if isempty(ah), ah=0;bh=0;ch=0; end 
%%%%%%%%% 
            [a,b,c,d]=series(ah,bh,ch,dh,a,b(:,1),-c(2,:),-d(2,1)); 
            pols=eig(a); 
            zers=tzero(a,b,c,d); 
% e1 -> y: HS  
          elseif chan==8 
%            [ah,bh,ch,dh,kh]=th2ss(th);         
%            [ah,bh,ch,dh]=minreal(ah,kh,ch,1);	corrected 21-10-1998: 
             [numh,denh]=th2tf(th,-1);[ah,bh,ch,dh]=tf2ss(numh,denh); 
             [ah,bh,ch,dh]=minreal(ah,bh,ch,dh); 
             if isempty(ah), ah=0;bh=0;ch=0; end 
%%%%%%%%% 
            [a,b,c,d]=series(ah,bh,ch,dh,a,b(:,2),c(2,:),d(2,2)); 
            pols=eig(a); 
            zers=tzero(a,b,c,d); 
          end 
          plot(real(pols),imag(pols),'*',real(zers),imag(zers),'o') 
          h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:)); 
          if ~isempty(zers),set(h(2),'color',modcol(mcact(i*j),:)),end 
          hold on 
        end 
      end 
    end 
% Get unit circle option 
    if strcmp(get(findobj(gcf,'label','Unit circle'),'checked'),'on')         
      om=linspace(0,2*pi);plot(cos(om),sin(om),'w') 
    end 
  
    set(fignum,'visible','on') 
  end 
%********************************************************************* 
% 
%              Input/output simulation 
% 
%********************************************************************* 
 
elseif option==3 
 % Only for parametric models 
 if strcmp(eval,'all') 
   fignum=findobj(0,'name',CLIDevas{option}); 
   figure(fignum) 
 % Calcultate input/output simulation 
    [DATA,x1,x2,p]=clobjget('datv'); 
    time=[0:1:length(DATA)-1]*p(2,1); 
    r2r1=DATA(:,[3 4]); 
  
   if strcmp(get(findobj(gcf,'label','Signal plot'),'checked'),'on')             
     subplot(211),cla,plot(time,DATA(:,1),'w-') 
     title('Measured  output (white) and simulated output y'),hold on 
     subplot(212),cla,plot(time,DATA(:,2),'w-') 
     title('Measured input (white) and simulated input u'), xlabel('time')   
     hold on 
   else 
     subplot(211),cla,title('Difference between measured and predicted output y'),hold on 
     subplot(212),cla,title('Difference between measured and predicted input u'),xlabel('time'),hold on 
   end  
for j=1:totcon 
  C=get(CLIDmods(21,conact(j)),'userdata'); 
    for i=1:totmod 
    % Only for parametric models 
    if CLIDmods(57,modact(i))==1 
      th=get(CLIDmods(11,modact(i)),'userdata'); 
%      [a,b,c,d]=th2ss(P); P=[a b;c d];  
      [numg,deng]=th2tf(th,1);[ag,bg,cg,dg]=tf2ss(numg,deng); 
      [ag,bg,cg,dg]=minreal(ag,bg,cg,dg); 
      if isempty(ag), ag=0;bg=0;cg=0; end 
      P=[ag bg;cg dg]; 
%%%%%%%%%%%%%%%%%%%%%%% 
      [TPC,nTPC]=cltpc(P,C,length(P)-1,length(C)-1); 
      [a,b,c,d]=clsplit(TPC,nTPC); 
      [a,b,c,d]=minreal(a,b,c,d);  
      YU=dlsim(a,b,c,d,r2r1); 
      if ~strcmp(get(findobj(gcf,'label','Signal plot'),'checked'),'on')             
       YU=DATA(:,[1 2])-YU; 
      end 
      subplot(211), plot(time,YU(:,1),'-') 
      h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:)); 
      hold on 
      subplot(212),plot(time,YU(:,2)) 
      h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:)); 
      hold on 
    end 
    end 
end 
 
    set(fignum,'visible','on') 
  % Get zoom option 
    if strcmp(get(findobj(gcf,'label','Zoom'),'checked'),'on')         
     zoom on 
    else 
     zoom off 
    end 
  end     
 
%********************************************************************* 
% 
%              Correlation test 
% 
%********************************************************************* 
 
elseif option==4 
 % Only for parametric models 
 if strcmp(eval,'all') 
   fignum=findobj(0,'name',CLIDevas{option}); 
   figure(fignum) 
 % Get number of lags for the correlation from user data from the figure  
   us=get(fignum,'userdata');N=us(2,1); 
   subplot(211),cla, 
   subplot(212),cla 
  
% Calculate correlation test 
for j=1:totcon 
    C=get(CLIDmods(21,conact(j)),'userdata'); 
    for i=1:totmod 
    % Only for parametric models 
    if CLIDmods(57,modact(i))==1 
      th=get(CLIDmods(11,modact(i)),'userdata'); 
      tsamp=gett(th);    %% Addition of 25/08/98 
%      [a,b,c,d]=th2ss(P); P=[a b;c d];  
      [numg,deng]=th2tf(th,1);[ag,bg,cg,dg]=tf2ss(numg,deng); 
      [ag,bg,cg,dg]=minreal(ag,bg,cg,dg); 
      if isempty(ag), ag=0;bg=0;cg=0; end 
      P=[ag bg;cg dg]; 
%%%%%%%%%%%%%%%%%%%%%%% 
      [TPC,nTPC]=cltpc(P,C,length(P)-1,length(C)-1); 
      [a,b,c,d]=clsplit(TPC,nTPC); 
      [a,b,c,d]=minreal(a,b,c,d);  
  %%Correction 25-08-98: 
	  DATA=clobjget('datv'); 
	  [num,den]=ss2tf(a,b(:,2),c(1,:),d(1,2)); 
      th=poly2th([],num,1,[],den,[],tsamp); 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%end correction		   
      [E,R]=clresid([DATA(:,1) DATA(:,4)],th,N); 
      subplot(211),hold on,h=get(gca,'children'); 
      set(h(1),'color',modcol(mcact(i*j),:)); 
      set(h(2),'color',modcol(mcact(i*j),:)); 
      set(h(3),'color',modcol(mcact(i*j),:)); 
  %%Correction 25-08-98: 
	  [num,den]=ss2tf(a,b(:,2),c(2,:),d(2,2)); 
      th=poly2th([],num,1,[],den,[],tsamp); 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%end correction		   
      [E1,R1]=clresid1([DATA(:,2) DATA(:,4)],th,N); 
      subplot(212),hold on,h=get(gca,'children'); 
      set(h(1),'color',modcol(mcact(i*j),:)); 
      set(h(2),'color',modcol(mcact(i*j),:)); 
      set(h(3),'color',modcol(mcact(i*j),:)); 
    end 
    end 
end 
 
    subplot(211) 
    zoom on,title('Cross corr. function between r and residual output y')  
    subplot(212) 
    zoom on,title('Cross corr. function between r and residual input u') 
    set(h,'visible','on') 
  end 
%********************************************************************* 
% 
%            Closed-loop step response 
% 
%********************************************************************* 
 
elseif option==5 
 % Only for parametric models 
 if strcmp(eval,'all') 
   fignum=findobj(0,'name',CLIDevas{option}); 
   figure(fignum) 
 % Get length of the step response and the channel from user data from the figure  
   us=get(fignum,'userdata'); 
   N=us(2,1);chan=us(1,1);                    
 
  % Get grid option 
    if strcmp(get(findobj(gcf,'label','Grid'),'checked'),'on')         
     xgrid='on';ygrid='on'; 
    else 
     xgrid='off';ygrid='off'; 
    end 
  % Set the figure layout  
   if chan==1|chan==2|chan==3|chan==4 
    subplot(1,1,1),cla,set(gca,'xgrid',xgrid,'ygrid',ygrid),ylabel('stepresponse'),xlabel('time'),hold on 
   elseif chan==5 
    subplot(221),cla,set(gca,'xgrid',xgrid,'ygrid',ygrid),hold on 
    subplot(222),cla,set(gca,'xgrid',xgrid,'ygrid',ygrid),hold on 
    subplot(223),cla,set(gca,'xgrid',xgrid,'ygrid',ygrid),xlabel('time'),hold on 
    subplot(224),cla,set(gca,'xgrid',xgrid,'ygrid',ygrid),xlabel('time'),hold on    
   end 
 
 % Calculate closed-loop step responses 
for j=1:totcon 
  C=get(CLIDmods(21,conact(j)),'userdata'); 
    for i=1:totmod 
    % Only for parametric models 
    if CLIDmods(57,modact(i))==1 
     % Calculate 
     th=get(CLIDmods(11,modact(i)),'userdata'); 
%      [a,b,c,d]=th2ss(th); P=[a b;c d];  
      [numg,deng]=th2tf(th,1);[ag,bg,cg,dg]=tf2ss(numg,deng); 
      [ag,bg,cg,dg]=minreal(ag,bg,cg,dg); 
      if isempty(ag), ag=0;bg=0;cg=0; end 
      P=[ag bg;cg dg]; 
%%%%%%%%%%%%%%%%%%%%%%% 
      [TPC,nTPC]=cltpc(P,C,length(P)-1,length(C)-1); 
      [a,b,c,d]=clsplit(TPC,nTPC); 
      [a,b,c,d]=minreal(a,b,c,d);  
     if strcmp(get(findobj(gcf,'label','Step response'),'checked'),'on')         
      stepresp1=dstep(a,b,c,d,1,N); 
      stepresp2=dstep(a,b,c,d,2,N);   
      ytext='step response'; 
     else 
      stepresp1=dimpulse(a,b,c,d,1,N); 
      stepresp2=dimpulse(a,b,c,d,2,N);   
      ytext='impulse response'; 
     end     
 % S: r1 -> u 
     if chan==1 
      subplot(1,1,1),plot(stepresp2(:,2)),title('S: r1 -> u'),xlabel('time'),ylabel(ytext) 
      h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:)); 
 % CS: r2 -> u 
     elseif chan==2 
      subplot(1,1,1),plot(stepresp1(:,2)),title('CS: r2 -> u'),xlabel('time'),ylabel(ytext) 
      h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:)); 
 % PS: r1 -> y 
     elseif chan==3 
      subplot(1,1,1),plot(stepresp2(:,1)),title('PS: r1 -> y'),xlabel('time'),ylabel(ytext) 
      h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:)); 
 % T: r2 -> y 
     elseif chan==4 
      subplot(1,1,1),plot(stepresp1(:,1)),title('T: r2 -> y'),xlabel('time'),ylabel(ytext) 
      h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:)); 
 % S: r1 -> u 
     elseif chan==5 
      subplot(221),plot(stepresp1(:,1)),title('T: r2 -> y'),ylabel(ytext) 
      h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:)); 
      subplot(222),plot(stepresp2(:,1)),title('GS: r1 -> y')  
      h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:));  
      subplot(223),plot(stepresp1(:,2)),title('CS: r2 -> u'),ylabel(ytext) 
      h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:)); 
      subplot(224),plot(stepresp2(:,2)),title('S: r1 -> u') 
      h=get(gca,'child');set(h(1),'color',modcol(mcact(i*j),:)); 
     end 
   end 
   end 
end 
    
    set(fignum,'visible','on') 
   % Get zoom option 
    if strcmp(get(findobj(gcf,'label','Zoom'),'checked'),'on')         
     zoom on 
    else 
     zoom off 
    end     
  end 
  
%********************************************************************* 
% 
%            Plant transfer function 
% 
%********************************************************************* 
 
 
elseif option==6 
 if strcmp(eval,'all') 
   fignum=findobj(0,'name',CLIDevas{option}); 
   figure(fignum) 
 
 % Get length of the step response and the channel from user data from the figure  
   us=get(fignum,'userdata');Nus=length(us); 
   wp=us(2:Nus,1);                   % Get frequency axis  
   chan=us(1,1); 
 
  % Get grid option 
    if strcmp(get(findobj(gcf,'label','Grid'),'checked'),'on')         
     xgrid='on';ygrid='on'; 
    else 
     xgrid='off';ygrid='off'; 
    end 
 
   % Clear the axis   
    subplot(211),cla,set(gca,'xgrid',xgrid,'ygrid',ygrid), 
    subplot(212),cla,set(gca,'xgrid',xgrid,'ygrid',ygrid), 
  
   % Calculate and frequency response of all selected models 
   for i=1:totmod 
    % Calculate  
     P=get(CLIDmods(11,modact(i)),'userdata'); 
     if CLIDmods(57,modact(i))==1 
      w=wp;  
      Gf=th2ff(P,1,w); 
      if chan==1 
       [w,m,f]=getff(Gf,1,1); 
      elseif chan==2 
        [w,m,f]=getff(Gf,0,1); 
        if isempty(m) 
          w=us(2:Nus,1);m=ones(size(w));f=zeros(size(w)); 
        end 
      end      
     elseif CLIDmods(57,modact(i))==0 
        [w,m,f]=getff(P); 
        if chan==2 
          m=ones(size(w));f=zeros(size(w)); 
        end         
     end     
% Get amplitude scale  
     if strcmp(get(findobj(gcf,'label','Log amplitude scale'),'checked'),'on')         
       yscale='log'; 
     else 
       yscale='lin'; 
     end 
 
     % Plot 
      subplot(211),loglog(w,m)  
      h=get(gca,'child');set(h(1),'color',modcol(modact(i),:)); 
	  %if yscale=='log', axis([min(w)/5 max(w)*5 1e-3 10]), end 
	  if chan==1 
	    title('Open-loop frequency response G: u -> y') 
      elseif chan==2 
	    title('Open-loop frequency response H: e1 -> y') 
	  end 
      hold on 
      subplot(212),semilogx(w,f) 
      h=get(gca,'child');set(h(1),'color',modcol(modact(i),:)); 
	  %if yscale=='log', axis([min(w)/5 max(w)*5 1e-3 10]), end 
	  hold on 
    end    
 
 
    set(fignum,'visible','on') 
    zoom on 
 end 
%********************************************************************* 
% 
%            Plant pole-zeros 
% 
%********************************************************************* 
 
elseif option==7 
  figure(fignum),cla 
% Get the channel from user data from the figure  
  us=get(fignum,'userdata');Nus=length(us); 
  chan=us(1,1); 
    
% Set figure layout right 
%  hx=axes('Position',[0.1 0.1 0.8 0.8]); % try-out.... 
 
% Get zoom option 
  if strcmp(get(findobj(gcf,'label','Zoom'),'checked'),'on')         
    zoom on 
  else 
    zoom off 
  end 
% Get grid option 
  if strcmp(get(findobj(gcf,'label','Grid'),'checked'),'on')         
    grid on 
  else 
    grid off 
  end 
% Get unit circle option 
  if strcmp(get(findobj(gcf,'label','Unit circle'),'checked'),'on')         
    om=linspace(0,2*pi);plot(cos(om),sin(om),'w') 
  end 
 
  if strcmp(eval,'all') 
    fignum=findobj(0,'name',CLIDevas{option}); 
% Calculate plant pole-zeros 
    for i=1:totmod 
% Only for parametric models 
      if CLIDmods(57,modact(i))==1 
        P=get(CLIDmods(11,modact(i)),'userdata'); 
        if chan==1 
          [zepo,k]=zp(P);   
        elseif chan==2 
          [zepo,k]=zp(P,0,1);   
        end 
        [ze,po]=getzp(zepo); 
        plot(real(ze),imag(ze),'o',real(po),imag(po),'*') 
        if [~isempty(ze) ~isempty(po)] 
          h=get(gca,'child');set(h([1 2]),'color',modcol(modact(i),:)); 
        end 
        hold on 
      end 
    end 
  end 
  set(fignum,'visible','on') 
end