function [AM,t,vx,vy,vw2,vz,vq,vxm,vqm]=simu(fz,fn,rz,rn,zs1,sdef,wv,zv)
%SIMU     step response of a single control loop - simulation and quality criteria
%         [AM,t,vx,vy,vw2,vz,vq,vxm,vqm]=simu(fz,fn,rz,rn,zs1,sdef,wv,zv)
%Outputs:        AM - analysis matrix=[hun t(hun) tar hue t(hue) taus;
%                                     RF1 RF2    RF3 RF4 RF5    0    ]
%  hun,hue,tar,taus - undershoot resp. overshoot; rise time resp. transient time
%  t(hun),t(hue)    - point in time of the occurrence of undershoot resp. overshoot
%  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 of error
%                     resp. IT2AE (quadratic time weighting)
%  t,vx,vy          - time, actual value, controller output
%  vw2,vz,vq        - filtered setpoint, disturbance variable, uncontrolled plant
%  vxm,vqm          - x resp. q mean values(MW), if mean value measuring
%Inputs:
%  fz,fn,rez,rn     - filter(f) and controller(r): numerator, denominator: z-domain
%  zs1              - polynomial between controller and disturbance action z-domain
%  sdef             - model of the continuous plant - optionally result of trans.m
%  wv,zv  - setpoint resp. disturbance variable vector, length=sampling_step_number
%Graphic representation: default colours depend on MATLAB version
%  cyan - setpoint(w2);     yellow - controller output; green - uncontrolled plant
%  blue - setpoint(20% wv); red - disturbance variable; magenta - actual value 
%  also valid for mean values

% KEY WORDS: simulation, quality criteria, single control loop

% Copyright (c) 1997-2023. All Rights Reserved.
%                     Dr. Gert-Helge Geitner, TU Dresden,
%                     https://wwwpub.zih.tu-dresden.de/~geitner/ae_1_0e.htm

[i1,Df]  = size(fz);  [i2,Bf] =size(fn);
[i3,D]   = size(rz);  [i4,B]  =size(rn);
[i5,i6]  = size(sdef); N2=sdef(1,1); M=sdef(2,1); MWM =sdef(3,4);
                       TA=sdef(3,1); m=sdef(3,2); kode=sdef(3,3);
[i7,ios1]= size(zs1); [i8,md] =size(wv);  [i9,j]=size(zv);
if ((i1+i2+i3+i4+i7+i8+i9)~=7)|(i5~=3)|...
   sum(abs(sdef(1,2+sdef(1,1):end)))~=0|sum(abs(sdef(2,2+sdef(2,1):end)))~=0
   error('Incorrect inputs'); end
if (md~=j) error('Input signal vectors of different length !'); end
for i=1:N2 zs2(i)=sdef(1,1+i); end; for i=1:M ns2(i)=sdef(2,1+i); end
iof=max(Df-1,Bf-1); ior=max(D-1,B-1); ios1=ios1-1;
ans=input('Output_step_size/sampling_time: standard(10:1)   yes(any)/no(1): ');
if isempty(ans)|ans~=1 Tk=10;
else
   Tk=input('Input Output_step_size/sampling_time= ');
   if ((Tk<1)|(rem(Tk,1)~=0)) error('Incorrect input'); end
end
gb=input('Graphic representation: standard=100% scaling   yes(any)/no(1): ');
nov=ones(1,4);
if isempty(gb)|gb~=1 gb=0;
else j=1;
   Te='Input weighting factors (x,y,w2,z) curve0 = ';
   while j<5
         Te(41)=48+j; Te=setstr(Te); i=input(Te);
         if      ((i<0)|(i>1)) disp('Incorrect input');
         elseif  i<0.1         disp('Not useful');
         elseif  isempty(i)    j=5;
         else nov(j)=i; j=j+1; end; end
end
if MWM
   MWAZ=input('Graphic representation of mean values x and q:  no(any)/yes(1): ');
   if isempty(MWAZ)|MWAZ~=1 MWAZ=0; end
else MWAZ=0; end
T=TA/Tk;j=fix((m*10^6)/(T*10^6)); m=rem(m*10^6,T*10^6)/10^6;
        % Die Erweiterung um "10^6" soll richtige Ergebnisse im wahr-
        % scheinlichen Zahlenbereich sichern, da "fix" und "rem" für
        % Zahlen < "1" getrennt betrachtet Fehler liefern können:
        % z.B.: fix(0.7/0.1)=6 und rem(0.7,0.1)=0.1
mdk=fix(j/Tk); Tkk=rem(j,Tk);  %Totzeiten: bez. Abtastzeit bzw. kont. Zeitauflösung
[z0,zn]=trans(zs2,ns2,T,m,kode-MWM);   N2=length(z0); % l=ceil(m/T), m=l*T-m
                                       M =length(zn); %M: totzeitunabhängig
if MWM==1
   [z0m,znm]=trans(zs2,ns2,TA,m,kode); N3=length(z0m);%im Normalfall N3=N2+1; aber:
                                       M3=length(znm); end%trans: Ordnungsgleichheit
j=0; while abs(z0(j+1))<(1e-14) j=j+1; end            %MWM: Totzeit j bleibt gleich
for i=(1+j):N2 zz(i-j)=z0(i); end               %für filter.m um Totzeit verschieben
    if MWM for i=(1+j):N3 zzm(i-j)=z0m(i); end; end
ios2=max(N2-j,M)-1; if MWM ios2m=max(N3-j,M3)-1; end  %Vektorlänge Anfangswerte
if iof ~=0 zif(iof)=0;   end
if ior ~=0 zir(ior)=0;   end
if ios1~=0 zis1(ios1)=0; end
if ios2~=0 zis2(ios2)=0;
           ziur(ios2)=0; end; x=0;
if MWM&ios2m~=0 zis2m(ios2m)=0;
                ziurm(ios2m)=0; end;
y2v(mdk+1)=0; sev(mdk+1)=0; vxv(Tkk+1)=0; vqv(Tkk+1)=0;
vx(1)=0;  vy(1)=0;  vw2(1)=0; vz(1)=0; vw(1)=0; vxw(1)=0; vq(1)=0;
if MWM vxm(1)=0; vqm(1)=0; xm=0; qm=0; end
for i=1:md
    if   iof~=0  [w2,zff] =filter(fz,fn,wv(i),zif);  zif=zff;
    else         [w2]     =filter(fz,fn,wv(i));      end
if MWM  xw=w2-xm;
else    xw=w2-x; end
    if ior~=0  [y,zfr]  =filter(rz,rn,xw,zir);       zir=zfr;
         else  [y]      =filter(rz,rn,xw);           end
    if ios1~=0 [y1,zfs1]=filter(zs1,1,y,zis1);       zis1=zfs1;
         else  [y1]     =filter(zs1,1,y);            end
    y2=y1-zv(i);    y2v(1)=y2; y2=y2v(1+mdk); %einschreiben,auslesen,schieben
    se=wv(i)-zv(i); sev(1)=se; se=sev(1+mdk);
       for j=1:mdk  y2v(mdk+2-j)=y2v(mdk+1-j); %Totzeit bezüglich Abtastzeit 
                    sev(mdk+2-j)=sev(mdk+1-j); end
    for j=1:Tk
        if ios2~=0 [x,zfs2]=filter(zz,zn,y2,zis2); zis2=zfs2;
                   [q,zfur]=filter(zz,zn,se,ziur); ziur=zfur;
             else  [x]     =filter(zz,zn,y2);
                   [q]     =filter(zz,zn,se);      end
        vxv(1)=x;     x=vxv(1+Tkk); %einschreiben,auslesen,schieben
        vqv(1)=q;     q=vqv(1+Tkk);
        for ii=1:Tkk  vxv(Tkk+2-ii)=vxv(Tkk+1-ii);     %Totzeit bezüglich konti-
                      vqv(Tkk+2-ii)=vqv(Tkk+1-ii); end %nuierlicher Zeitauflösung
        vx(2+j+(i-1)*(Tk+1)) =x;
        vy(1+j+(i-1)*(Tk+1)) =y;
        vw2(1+j+(i-1)*(Tk+1))=w2;
        vz(1+j+(i-1)*(Tk+1)) =zv(i);
        vw(1+j+(i-1)*(Tk+1)) =wv(i);
        vxw(1+j+(i-1)*(Tk+1))=xw;
        vq(2+j+(i-1)*(Tk+1)) =q;
        if MWM vxm(1+j+(i-1)*(Tk+1))=xm;
               vqm(1+j+(i-1)*(Tk+1))=qm; end
    end
    if MWM
      if ios2m~=0 [xm,zfs2m]=filter(zzm,znm,y2,zis2m); zis2m=zfs2m;
                  [qm,zfurm]=filter(zzm,znm,se,ziurm); ziurm=zfurm;
      else        [xm]      =filter(zzm,znm,y2);
                  [qm]      =filter(zzm,znm,se);      end
    end
    vx(3+j+(i-1)*(Tk+1)) = vx(2+j+(i-1)*(Tk+1));% Treppenfunktn ermöglichen
    vy(2+j+(i-1)*(Tk+1)) = vy(1+j+(i-1)*(Tk+1));%       "            "
    vw2(2+j+(i-1)*(Tk+1))=vw2(1+j+(i-1)*(Tk+1));%       "            "
    vz(2+j+(i-1)*(Tk+1)) = vz(1+j+(i-1)*(Tk+1));%       "            "
    vw(2+j+(i-1)*(Tk+1)) = vw(1+j+(i-1)*(Tk+1));%       "            "
    vxw(2+j+(i-1)*(Tk+1))=vxw(1+j+(i-1)*(Tk+1));%       "            "
    vq(3+j+(i-1)*(Tk+1)) = vq(2+j+(i-1)*(Tk+1));%       "            "
    if MWM
       vxm(2+j+(i-1)*(Tk+1))=vxm(1+j+(i-1)*(Tk+1));%    "            "
       vqm(2+j+(i-1)*(Tk+1))=vqm(1+j+(i-1)*(Tk+1));end% "            "
end
 vw2(3+j+(i-1)*(Tk+1))=vw2(2+j+(i-1)*(Tk+1));% Zeitvektor auffüllen
 vz(3+j+(i-1)*(Tk+1)) = vz(2+j+(i-1)*(Tk+1));%     "          "
 vw(3+j+(i-1)*(Tk+1)) = vw(2+j+(i-1)*(Tk+1));%     "          "
 if MWM  xw=w2-xm;
 else    xw=w2-x;  % w2-Brechnung nicht nötig, mit Filter keine Auswertung
 end
 if ior~=0  [y,zfr]  =filter(rz,rn,xw,zir);
    else    [y]      =filter(rz,rn,xw); end
 vxw(3+j+(i-1)*(Tk+1))=xw;
  vy(3+j+(i-1)*(Tk+1))=y;
 if MWM vxm(3+j+(i-1)*(Tk+1))=xm;
        vqm(3+j+(i-1)*(Tk+1))=qm; end
szma=max([max(wv) max(zv)]); szmi=min([min(wv) min(zv)]);
if szma<0 szma=0; end
if szmi>0 szmi=0; end
wabs=max([szma abs(szmi)]); vwn = vw/wabs*100;
                            vw2n=vw2/wabs*100;
                            vxn = vx/wabs*100;
                            vzn = vz/wabs*100;
                           i=2; while sdef(1,i)==0; i=i+1; end; Vp=sdef(1,i);
                            vqn = vq/(wabs*Vp)*100;
                            if MWM vqmn=vqm/(wabs*Vp)*100;
                                   vxmn=vxm/wabs*100; end
snmi=szmi/wabs*100; snma=szma/wabs*100;
wema=max([max(vw2n) max(vxn)]); wemi=min([min(vw2n) min(vxn)]);
if (snma-wema)<0 yma=snma-(snma-wema); else yma=snma; end
if (snmi-wemi)>0 ymi=snmi-(snmi-wemi); else ymi=snmi; end

stma=max(vy); if stma<0 stma=0; end
stmi=min(vy); if stmi>0 stmi=0; end
wabs=max([stma abs(stmi)]); vyn = vy/wabs*100;
stnma=stma/wabs*100; stnmi=stmi/wabs*100;
if stnma>yma yma=stnma; end
if stnmi<ymi ymi=stnmi; end
if gb vxn=vxn*nov(1); vyn=vyn*nov(2); vw2n=vw2n*nov(3); vzn=vzn*nov(4);
   if MWM vxmn=vxmn*nov(1); end; end
axis([0 md ymi*1.05 yma*1.05])
t=-1:Tk-1; t(1)=0;           % Darstellung von Treppenkurven ermöglichen
for i=1:md-1 for j=1:Tk+1 t((Tk+1)*i+j)=t(j)+Tk*i; end, end
t(1+md*(Tk+1))=t(md*(Tk+1))+1; t(2+md*(Tk+1))=t(1+md*(Tk+1)); t=t/Tk;
%4.2c:
%plot(t,vxn,'c11',t,vyn,'c5',t,vw2n,'c12',t,vzn,'c2',t,vwn*0.2,'c10',...
%     t,vqn,'c3');
%höher:
if ~MWAZ plot(t,vqn,'g',t,vyn,'y',t,vw2n,'c',t,vzn,'r',t,vwn*0.2,'b',...
              t,vxn,'m');
else     plot(t,vqn,'g',t,vyn,'y',t,vw2n,'c',t,vzn,'r',t,vwn*0.2,'b',...
              t,vxn,'m',t,vxmn,'m',t,vqmn,'g');
end
xlabel('time/samplig_time'); ylabel('per_cent'); grid;
title('Please note: representation is normalized!');
disp('Stop: please press any key.'); pause
vxR=vx; vyR=vy;           %Retten wegen nachfolgender Abfrage zum Ausgangssetzen
if MWM vxmR=vxm; end
ans=input('Outputs = simulation_data(SD) or drawings(DR) ? SD(any)/DR(1): ');
if isempty(ans) ans=0; end
if ans==1 vx=vxn; vy=vyn; vw2=vw2n; vz=vzn; vq=vqn;
   if MWM vxm=vxmn; vqm=vqmn; end; end
ans=input('Analyses of the response behaviour no(any)/yes(1): ');
if isempty(ans) ans=0; end
if ans~=1 AM=0; if ~MWM vxm=0; vqm=0; end
else TvP=1;  %Position des letzten Zeitwertes vor dem Sprung
   while ((wv(TvP)==0)&(zv(TvP)==0)&(TvP<md)) TvP=TvP+1; end; bte=0;
   if zv(TvP)~=0 bfu=0; prwe=zv(TvP);
                 if sum(abs(wv))~=0|sum(abs(zv))~=abs(sum(zv)); bte=1; end
   else          bfu=1; prwe=wv(TvP);
                 if sum(abs(zv))~=0|sum(abs(wv))~=abs(sum(wv)); bte=1; end
   end
  Tv=t(1+(Tk+1)*(TvP-1)); %normierter Zeitwert bei Sprung
  if bte disp([07 07 'Prepared for a single stepped response only !']);
  else
   vxn=vxn/nov(1); vyn=vyn/nov(2); vw2n=vw2n/nov(3); vzn=vzn/nov(4);
   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('Input range of tolerance (in %) = ');
            if  toba>10|toba<0|isempty(toba) disp('Illegal value');
            else i=1; end; end
   end;
   disp('Choice of a sampling integral of error - (0=all,Enter=no): ');
  brf=input('quadratic-1,absolute_linear-2,enhanced quadratic-3,ITAE-4,IT2AE-5: ');
   if ((brf<0)|(brf>5)|(rem(brf,1)~=0)) error('Incorrect input'); end
   if isempty(brf) brf=6; end
   if ((brf==5)|(brf==0))
      ans=input('Weighted controller output: standard=0.1   yes(any)/no(1): ');
      if isempty(ans)|ans~=1 lampda=0.1;
      else lampda=input('Input - controller output weighting = ');
         if lampda<0|lampda>1|isempty(lampda) error('Incorrect input'); end;
      end; end
   if MWM
        MWAW=input('Quality criteria by instantaneous_value(any)/mean_value(1): ');
        if isempty(MWAW)|MWAW~=1 MWAW=0; else vxn=vxmn/nov(1); end
   else MWAW=0; end
   disp('Quality criteria of the selected calculation intreval result in:');
   if prwe<0 vxn=-vxn; vwn=-vwn; vzn=-vzn; end
   [hun,t1P]=min(vxn); t1=t(t1P)-Tv;
   if abs(hun)<0.01 hun=0; t1=0; disp('Undershoot represents 0 %');
   else disp('Undershoot in unit % at time_t/sampling_time:');
    disp([hun t1]); end
   [hue,t2P]=max(vxn); t2=t(t2P)-Tv;
   if bfu hue=hue-100; vgv=vwn; pep=1+toba/100; mep=1-toba/100;
   else                vgv=vzn; pep=toba/100;   mep=-pep; end
   if hue<0.01 hue=0; t2=0; disp('Overshoot represents 0 %');
   else disp('Overshoot  in unit % at time_t/sampling_time:');
     disp([hue t2]); end
   i=t1P; an=2+md*(Tk+1); j=an+1; %Störsprung: eintauchen statt verlassen erkennen
   while i<an
         if  (vxn(i+1)<=(vgv(i+1)*mep)) i=i+1; else j=i+1; i=an; end; end
   if j==an+1 tar=inf; else tar=t(j)-Tv; end
   disp('Rise_time/sampling_time = '); disp(tar);
   i=0; j=an;
   while i<an
         if ((vxn(an-i)<=(vgv(an-i)*pep))&(vxn(an-i)>=(vgv(an-i)*mep)))
         i=i+1; else j=an-i+1; i=an; end; end
   if    j==an+1; taus=inf; elseif j==an taus=0; else taus=t(j)-Tv; end
   disp('Transient_time/sampling_time = '); disp(taus);
   AM(1)=hun; AM(2)=t1; AM(3)=tar; AM(4)=hue; AM(5)=t2; AM(6)=taus;
   if brf~=6
   RF1=0; RF2=0; RF3=0; RF4=0; RF5=0; vyn=vyR/prwe;
   if MWAW vxw=(vw-vxmR)/prwe; else vxw=(vw-vxR)/prwe; end
   tn=t; k=0;                                                  %Doppelung streichen
   for i=1:Tk+1:an tn(i-k)=[]; vxw(i-k)=[]; vyn(i-k)=[]; k=k+1; end;
   Anz=an-k;                  %verringerte Werteanzahl wegen gestrichener Doppelung
   Bg=1+(TvP-1)*Tk; %Korrektur für Berechnungsbeginn bei verzögertem Eingangssignal
if (brf==2)|(brf==0)
   for i=Bg:Anz RF2=RF2+abs(vxw(i)); end; end
if (brf==1)|(brf==0)
   for i=Bg:Anz RF1=RF1+(vxw(i))^2; end; end
if (brf==4)|(brf==0)
   for i=Bg:Anz RF4=RF4+(tn(i)-Tv)*abs(vxw(i)); end; RF4=RF4*TA; end %Echtzeitbezug
if (brf==5)|(brf==0)
   for i=Bg:Anz RF5=RF5+(tn(i)-Tv)^2*abs(vxw(i)); end; RF5=RF5*TA^2; end %   "
if (brf==3)|(brf==0)
   for i=Bg:Anz RF3=RF3+((vxw(i))^2)+lampda*((vyn(i)-vyn(Anz))^2); end; end;
   vRF=[RF1 RF2 RF3 RF4 RF5]*(TA/Tk); %konstante lineare Delta-t-Auflösung beachten
   if brf==0 disp('Integrals of error represent:');
      format short e; disp(vRF); format; brf=5;
   else disp('The selected integral of error represents:'); disp(vRF(brf)); end
   for i=1:brf AM(2,i)=vRF(i); end
   end
  end
end
