function [AM,tn,vwn,vxn]=guete(T,t,vw,vx,vy,vz,SV,UVA)
%GUETE    calculation of step response quality criteria using SIMULINK simulation
%         results of a closed loop
%         [AM(,tn,vwn,vxn)]=guete(T,t,vw,vx,vy,vz(,SV,UVA))
% Inputs:
%   T     - sampling time
%   t     - time                                    not normalized (column vector)
%   vw,vx - setpoint          resp. actual value      "      "        "      "
%   vy,vz - controller output resp. disturbance var.  "      "        "      "
%   (,SV          - optionally menu or control vector driven; control vector: 
%                   SV(1): range of tolerance in per cent (0<= valid tol. <=10) 
%                   SV(2): 1=RF1, 2=RF2, 3=RF3, 4=RF4, 5=RF5; 0=all, 6=none
%                          default=6 (active only for setting SV(1)!)
%                   SV(3): controller output weighting (COW), (0<= COW <=1),
%                          required only if SV(2)=3|0, default=0.1
%    UVA)         - optionally elimination of text outputs, e.g. in case one-on-one
%                   result outputs, or in case of using menu (SV=[]): any=no/1=yes
% Outputs:
%   AM            - analysis matrix=[hun t(hun) tar hue t(hue) taus;
%                                    RF1 RF2    RF3 RF4 RF5    0    ]
%   hun,hue       - undershoot resp. overshoot
%   t(hun),t(hue) - point in time of the occurrence of undershoot resp. overshoot
%   tar,taus      - rise time resp. transient time
%   RF1,RF2       - quadratic resp. absolute linear sampling integral of error
%   RF3           - enhanced quadratic sampling integral of error
%                   (additional weighted square of the controller output)
%   RF4,RF5       - ITAE (Integral of Time And Error) - sampling integral
%                   resp. IT2AE (quadratic time weighting) - sampling integral
%   (,tn,vwn,vxn) - optionally outputs: normalized vectors t,vw,vx

%  KEY WORDS: quality criteria, rise time, transient time, overshoot, undershoot,
%             integral of error, range of tolerance, SIMULINK

%  Copyright (c) 1997-2018. All Rights Reserved. 
%                     Dr. Gert-Helge Geitner; TU Dresden, Fak. ET,
%                     Elektrotechnisches Institut (ETI); Mommsenstr. 13;
%                     D-01062 Dresden, Germany;
%                     http://eeiwzg.et.tu-dresden.de/ae2_files/ae_1.htm

%---------------------------German-------------------------------------------------
%GUETE    Gütekenngrößen der Übergangsfunktion nach Regelkreissimulation (SIMULINK)
%         [AM(,tn,vwn,vxn)]=guete(T,t,vw,vx,vy,vz(,SV,UVA))
% Eingänge:
%   T       - Abtastzeit
%   t       - Zeit                  nicht normiert (Spaltenvektor)
%   vw,vx   - Soll-  bzw. Istwert     "      "           "
%   vy,vz   - Stell- bzw. Störgröße   "      "           "
%   (,SV          - wahlweise Abschaltung der Menüführung; Steuervektor: 
%                   SV(1): Toleranzband in % (0<= zul.Tol <=10) 
%                   SV(2): 1=RF1, 2=RF2, 3=RF3, 4=RF4, 5=RF5; 0=alle, 6=keine
%                          default=6 (nur wirksam, wenn SV(1) vorgegeben!)
%                   SV(3): Stellgrößenwichtung, falls SV(2)=3|0; (0<= SGW <=1)
%                          default=0.1
%    UVA)         - wahlweise Unterdrückung von Ausgabetexten, z.B. Einzel-  
%                   ergebnisausgabe (auch mit Menü bei SV=[]): bel.=nein/1=ja
% Ausgänge:
%   AM      - Auswertematrix=[hun t(hun) tar hue t(hue) taus;
%                             RF1 RF2    RF3 RF4 RF5    0    ]
%   hun,hue       - Unter- bzw. Überschwingweite
%   t(hun),t(hue) - Zeitpunkt des Auftretens von hun bzw.hue
%   tar,taus      - An- bzw. Ausregelzeit
%   RF1,RF2       - Quadratische bzw. Betragslineare Abtast-Regelfläche
%   RF3           - Erweitert quadratische Abtast-Regelfläche
%                   (gewichtetes Stellgrößenquadrat)
%   RF4,RF5       - ITAE bzw. IT2AE (quadratische Zeitwichtung) Abtast-RF
%   (,tn,vwn,vxn) - wahlweise Ausgänge: normierte Vektoren t,vw,vx

%  SCHLÜSSELWORTE: Regelgüte, Anregelzeit, Ausregelzeit, Überschwingweite,
%                  Integralkriterien, SIMULINK, Toleranzband

%  Copyright (c) 1997-2018. All Rights Reserved. 
%                     Dr. Gert-Helge Geitner; TU Dresden, Fak. ET,
%                     Elektrotechnisches Institut (ETI); Mommsenstr. 13;
%                     D-01062 Dresden, Germany;
%                     http://eeiwzg.et.tu-dresden.de/ae2_files/ae_1.htm
%-----------------------End German----------------------------------------------

TBMI=0; TBMA=10;    %range of tolerance: limits
LAMI=0; LAMA=1;     %Lampda: limits
FT01='Error: control vector SV incorrect';
FT02='Error: inputs incorrect';
FT03='Error: function "Guete" prepared for one step response only!';
KT01='Selected time interval - results:';
KT02='Undershoot represents 0%';
KT03='Undershoot in unit % at time_t/sampling_time:';
KT04='Overshoot represents 0%';
KT05='Overshoot  in unit % at time_t/sampling_time:';
KT06='Rise_time/sampling_time = ';
KT07='Transient_time/sampling_time = ';
KT08='Integrals of error represent:';
KT09='The selected integral of error represents:';

narginchk(6,8);   %obsolete versions: error(nargchk(6,8,nargin));
if ~(nargin==8) UVA=0; else if isempty(UVA)|~UVA==1; UVA=0; end; end
if nargin==8&isempty(SV)  MUEF=1; end
if nargin>6&~isempty(SV) [i1,j1]=size(SV); MUEF=0;
   if i1~=1|j1>3 error(FT01); end;
   toba=SV(1); if ((toba>TBMA)|(toba<TBMI)) error(FT01); end
   if j1==1 brf=6; else brf=SV(2);
                        if ((brf<0)|(brf>6)|(rem(brf,1)~=0)) error(FT01); end;
   end
   if j1==2&(brf==0|brf==3) lampda=0.1;
   elseif j1==3 lampda=SV(3); if ((lampda<LAMI)|(lampda>LAMA)) error(FT01); end;
   end
else MUEF=1; end                                                  
t=t'; vw=vw'; vx=vx'; vy=vy'; vz=vz';
[i,j]=size(T);    [i1,j1]=size(t);   [i2,j2]=size(vw);
[i3,j3]=size(vx); [i4,Anz]=size(vy); [i5,j5]=size(vz);
if i~=1|i1~=1|i2~=1|i3~=1|i4~=1|i5~=1|j~=1|T<=0|...
                                         ((j1+j2+j3+Anz)/4~=j5) error(FT02); end
%Im Vergleich zu simu.m: keine normierte Zeit, Tk nicht existent, T als Eingabe
TvP=1;                          %Position des letzten Zeitwertes vor dem Sprung
while ((vw(TvP)==0)&(vz(TvP)==0)&(TvP<Anz)) TvP=TvP+1; end; bte=0;
if vz(TvP)~=0 bfu=0; prwe=vz(TvP);
              if sum(abs(vw))~=0|sum(abs(vz))~=abs(sum(vz)); bte=1; end
else          bfu=1; prwe=vw(TvP);
              if sum(abs(vz))~=0|sum(abs(vw))~=abs(sum(vw)); bte=1; end; end;
Tv=t(TvP); %Zeitwert bei Sprung
if bte disp([07 07 FT03]); AM=0;t=0;vwn=0;vxn=0; return %Disp-Fehleranzeige
else
 if MUEF==1
      ans=input('range of tolerance: standard=1%   yes(any)/no(1): ');
      if isempty(ans)|ans~=1 toba=1;
      else i=0;
           while i==0 toba=input('range of tolerance (in %) = ');
                 if  ((toba>TBMA)|(toba<TBMI)) disp('Illegal value');
                 else i=1; end; end; end;
   disp('Choice of a sampling integral of error - (0=all,Enter=no): ');
   i=0; while i==0
  brf=input('Quadratic-1,absolute_linear-2,enhanced_quadratic-3,ITAE-4,IT2AE-5: ');
   if isempty(brf) brf=6; end
   if ((brf<0)|(brf>6)|(rem(brf,1)~=0)) disp('Incorrect input'); else i=1; end
   end
   if ((brf==3)|(brf==0))
      ans=input('Controller output weighting: standard=0.1   yes(any)/n0(1): ');
      if isempty(ans)|ans~=1 lampda=0.1;
      else i=0;
           while i==0 lampda=input('Input of controller output weighting = ');
              if ((lampda<LAMI)|(lampda>LAMA)) disp('Incorrect input');
              else i=1; end; end; end; end
 end
   if UVA==0 disp(KT01); end
   vxn=vx*100/prwe; vwn=vw*100/prwe; %vwn
   [hun,t1P]=min(vxn); t1=(t(t1P)-Tv); t1=t1/T; %Normierung
       if abs(hun)<0.01 hun=0; t1=0; if UVA==0 disp(KT02); end
       else if UVA==0 disp(KT03); disp([hun t1]); end; end
   [hue,t2P]=max(vxn); t2=(t(t2P)-Tv); t2=t2/T; %Normierung
       if bfu hue=hue-100; vgv=vwn; pep=1+toba/100; mep=1-toba/100;
       else                vgv=vz*100/prwe; pep=toba/100;   mep=-pep; end
       if hue<0.01 hue=0; t2=0; if UVA==0 disp(KT04); end
       else if UVA==0 disp(KT05); disp([hue t2]); end; end
   i=t1P; j=Anz+1; %t1P: bei Störsprung eintauchen statt verlassen von TB erkennen
   while i<Anz
         if  (vxn(i+1)<=(vgv(i+1)*mep)) i=i+1; else j=i+1; i=Anz; end; end
   if j==Anz+1 tar=inf; else tar=(t(j)-Tv)/T; end
      if UVA==0 disp(KT06); disp(tar); end
   i=0; j=Anz;
   while i<Anz
         if ((vxn(Anz-i)<=(vgv(Anz-i)*pep))&(vxn(Anz-i)>=(vgv(Anz-i)*mep)))
         i=i+1; else j=Anz-i+1; i=Anz; end; end
   if j==Anz+1; taus=inf; elseif j==Anz taus=0; else taus=(t(j)-Tv)/T; end
      if UVA==0 disp(KT07); disp(taus); end
   AM(1)=hun; AM(2)=t1; AM(3)=tar; AM(4)=hue; AM(5)=t2; AM(6)=taus;
%Wegen möglicher variabler Schrittweite sind im Vergleich zu simu.m keine 
%äquidistanten Zeitschritte voraussetzbar:
   if brf~=6
   RF1=0; RF2=0; RF3=0; RF4=0; RF5=0; vy=vy/prwe; vxw=(vw-vx)/prwe;
   for i=TvP:Anz-1 delta=t(i+1)-t(i);
    if (brf==2)|(brf==0) RF2=RF2+delta*abs(vxw(i)); end
    if (brf==1)|(brf==0) RF1=RF1+delta*(vxw(i))^2; end
    if (brf==4)|(brf==0) RF4=RF4+delta*(t(i)-Tv)*abs(vxw(i)); end
    if (brf==5)|(brf==0) RF5=RF5+delta*(t(i)-Tv)^2*abs(vxw(i)); end
   if (brf==3)|(brf==0) RF3=RF3+delta*(((vxw(i))^2)+lampda*((vy(i)-vy(Anz))^2));end
   end   
   sRF=abs('RF0');
   for i=1:5 AM(2,i)=eval(char(sRF+[0 0 i])); end
   if brf==0 if UVA==0 disp(KT08); format short e; disp(AM(2,1:5)); format; end;
   else if UVA==0 disp(KT09); disp(AM(2,brf)); end; end
   end
  end
vxn=vxn*sign(prwe); vwn=vwn*sign(prwe); tn=t/T;
