function [out1,out2,out3,out4] = eez(inp1,inp2,inp3,varargin)
%EEZ      Calculation of a discontinuous controller using the method of
%         Finite Time Settling (EEZ=EndlicheEinstellZeit) - starting via a plant
%         given in z-domain, typically applying function "trans" in toolbox BOD.
%         LTI syntax optionally; menu driven optionally;
% ==> There are two input versions valid:
%     a, standard MATLAB input syntax: numerator and denominator polynomial
%        [out1,out2(,out3,out4)] = eez(inp1,inp2(,inp3,varargin))
%                     inp1 - numerator   of the z-transformed plant
%                     inp2 - denominator of the z-transformed plant
%        Leading zeros of "inp1" and zeros at the end of "inp2" are evaluated!
%        Denominator(1)=1 is expected!
%     b, LTI model MATLAB syntax:
%        [out1(,out2)] = eez(inp1(,inp2,[],varargin))
%                        inp1 - LTI model of the z-transformed SISO plant
% ==> According input syntax there are two versions for the output of results:
%     a, standard syntax:
%             out1 - controller numerator   polynomial (descending powers of z)
%             out2 - controller denominator polynomial (    "         "   "  ")
%             out3 - pre-filter numerator   polynomial (    "         "   "  ")
%             out4 - pre-filter denominator polynomial (    "         "   "  ")
%     b, LTI syntax: out1 - LTI model - controller
%                    out2 - LTI model - pre-filter
% ==> (,out3,out4) resp. (,out2) - for delayed input signals useful only.
% ==> inpx = inp3(standard syntax) resp. inp2(LTI syntax) = menu switch off via
%                                                       following control
%                                                       vector:
% ==> inpx = inp3(standard syntax) resp. inp2(LTI syntax) = menu switch off via
%                                                     following control vector:
%     inpx(1) - any/1 = Optimization regarding: undelayed / delayed input signals
%     inpx(2) - -1/number = pole compensation if undelayed inputs with pre-filter
%               yes all poles(-1) / compensate a limited number (>0) / no(0)
%     inpx(3) - -1/number = zeros compensation if delayed input signals:
%               yes all zeros(-1)/ compensate a limited number (>0) / no(0)
%     inpx(4) - 0/number = undelayed input: additional/suboptimal rise time steps
%               0: none / >0: automatic- / <0: input- setting of coefficients
%     inpx(5) - 0/number = delayed input: additional/suboptimal rise time steps
%               0: none / >0: automatic- / <0: input- setting of coefficients
%     inpx(6) - any/1/2 = elimination of text outputs: no / yes(C) / yes(W+C)
%               (W - warnings; C - comments)
%     inpx(7) - any/1 = additional Difference Equation Element (DEE): no/yes
%               undelayed inp.: insert between controller and continuous plant
%               delayed inp.: insert between controller and disturbance action
%     inpx(8) - 0/1/2 = consideration of the input signal wave form by
%               wave form(0) / number integral parts(1) / setpoint polynomial(2)
%     inpx(9) - 0/1/2           if inpx(8)=0: step(0) / ramp(1) / parabola(2)
%               order           if inpx(8)=1: order of the integral part
%               numerator order if inpx(8)=2: input signal numerator polynomial 
%     inpx(10)- denominator order if inpx(8)=2: input signal denominator polynom.
%     inpx(11)- 0/number = consideration of integral parts of the plant to the
%               controller design: no(0)/yes(number of plants integral parts)
%
% ==> varargin useful only, if ~isempty(inpx) >>>AND<<< following conditions
%     AFu vector - if inpx(4)<0: external polynomial coefficients (descending) 
%                  undelayed inputs: suboptimal rise time steps: [afu1 afu2 ...]
%     ASt vector - if inpx(5)<0: external polynomial coefficients (descending) 
%                  delayed inp.: suboptimal transient time steps: [ast1 ast2 ...]
%     DGG vector - if inpx(7)=1: parameters - Difference Equation Element (DEE):
%                  [gain dead_time order coefficient1 coefficient2 ...]
%     SWZ vector - if inpx(8)=2: external polynomial coefficients (descending) 
%                  regarding setpoint numerator polynomial: [swz1 swz2 ...]
%     SWN vector - if inpx(8)=2: external polynomial coefficients (descending) 
%                  regarding setpoint denominator polynomial: [swn1 swn2 ...]
%     >>>ATTENTION<<< please do follow the sequence and do use empty vectors if 
%                     necessary.
% any == any other input, e.g. enter or 0, matches preferential input setting;
% >>The elimination of text outputs is recommended after getting familiar only!<<
% Defaults and completions of missing elements of inpx: zero

% KEY WORDS: Finite Time Settling control, discontinuous control, control loop,
%            optimization regarding undelayed input signals, optimization 
%            regarding delayed input signals, pre-filter design, optionally pole
%            compensation, optionally menu driven, optionally LTI syntax,
%            suboptimal solution, equalizing controller outputs, input signal
%            order

%  Copyright (c) 1997-2012. 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
%  Reference regarding the application of Finite Time Settling: Geitner, G.-H.:
%                     Entwurf digitaler Regler für elektrische Antriebe.
%                     publisher: VDE-Verlag GmbH Berlin und Offenbach 1996
%                     ISBN 3-8007-1847-2

FT01='Incorrect input';
FT02='No useful input';
FT03='Application to systems with feed through illegal!';
FT04='This plant features no zeros ~=0';
FT05='This plant features no poles';
FT06='Input argument error';
FT07='Input argument control_vector(inpx) - incorrect dimension';
FT08='Input argument control_vector(inpx) - element 4 or/and 5 incorrect!';
FT09='Incorrect DGG vector(varargin3) or disagreement with inpx(7)';
FT10='Incorrect ASt vector(varargin2) or disagreement with inpx(5)';
FT11='Incorrect AFu vector(varargin1) or disagreement with inpx(4)';
FT12='Deliver the coefficients of the plant as vectors!';
FT13='Denominator_of_the_plant(1)=1 expected!';
FT14='Incorrect SWZ vector(varargin4) or disagreement with inpx(9)';
FT15='Incorrect SWN vector(varargin5) or disagreement with inpx(10)';
FT16='Disagreement with inpx(8) and SWZ resp. SWN vector(varargin4/5)!';
FT17='Input argument control_vector(inpx) - element 8 incorrect!';
FT18='Input argument control_vector(inpx) - element 9 incorrect!';
FT19='Input argument control_vector(inpx) - element 11 incorrect!';
FT20='Control_vector(inpx) - element 5 incorrect or not useful!';
FT21='Control_vector(inpx) - element 4 incorrect or not useful!';
WT01='Warning: order of input vectors > 9 in fact?';
WT02='Warning: dead time > 3 in fact?';
WT03='Warning: the result matches an approximate solution!!!';
WT04='Warning: automatic increase of the number of compensated zeros!';
WT05='Warning: minimum controller order via full Nst compensation only!';
WT06='Warning: automatic increase of the number of compensated poles!';
WT07='Warning: minimum controller order via full pole compensation only';
WT08='Warning: minimum pre-filter order via full pole compensation only';
WT09='Warning: dead_time=T and factor T are considered as factored out!';
WT10='Warning: dead_time=T and factor T^2/2 are considered as factored out!';
KT01='Precondition of delayed input signals: mixing point at the plants input';
KT02='Consideration of the input signal waveform by choice of:';
KT03='Number of the integral parts of the plant Is=';
KT04='Number of the integral parts to be taken into account:';
KT05='Definition of a finite difference equation: Z1=Vz1*(1+g1*z^-1+...)*z^-k1';
KT06='The plant features following zeros:';
KT07='Compensation of the zero';
KT08='Minimum number of rise time steps for delayed input signals: ';
KT09='Automatic calculation of Ast polynomial via equalizing controller outputs';
KT10='The plant features following poles:';
KT11='Compensation of';
KT12='Minimum number of rise time steps for undelayed input signals: ';
KT13='Automatic calculation of Afu polynomial via equalizing controller outputs';
IT01='Delayed input optim. by zero compensation: yes(any)/no(1)=';
IT02='Delayed input: number additional/suboptim. rise time steps (Enter=0) =';
IT03='Setting of As polynomial coefficients: automatic(any)/input(1):';
IT04='Undelayed input optim. by pole compensat.: yes(any)/no(1)';
IT05='Prefilter design by pole compensat.: yes(any)/no(1) =';
IT06='Undelayed Inp.: number addition./suboptim. rise time steps (Enter=0)=';
IT07='Setting of Af polynomial coefficients: automatic(any)/input(1):';

if (nargin<1)|(nargin>8)|nargin==7|(nargin==2&isa(inp2,'lti'))|...
   (nargin==3&isa(inp3,'lti')) error(FT06); end
if isa(inp1,'lti')          %LTI-Modell übergeben
  if (nargin==3)|(nargin>3&~isempty(inp3))|~issiso(inp1)|~isdt(inp1) error(FT06);
   end
   if ~isa(inp1,'tf') inp1=tf(inp1); end
   [n,ns2]=tfdata(inp1,'v'); LTI=1;
   if nargin>1 MUEF=0; STEUV=inp2; else MUEF=1; STEUV(6)=0; end
else                        %alte Syntax übergeben
   if (nargin<2) error(FT06); end
   n=inp1; ns2=inp2; LTI=0;
   if nargin>2 MUEF=0; STEUV=inp3; else MUEF=1; STEUV(6)=0; end
end
[znum,snum]=size(STEUV); if znum>1|snum>11 error(FT07); end
if snum<11 STEUV(11)=0; end %Default-Ergänzung
if rem(STEUV(4),1)~=0|rem(STEUV(5),1)~=0 error(FT08); end
           %für restliche STEUV-Prüfung kann Eingabeprüfung verwendet werden!
if nargin>=6 [afu,ast,DGG]=deal(varargin{1:3});
   if nargin==8 [zqw,nqw]=deal(varargin{4:5}); [i,j]=size(zqw);[i1,j1]=size(nqw);
      if STEUV(8)==2&(i==0|i1==0) error(FT16); end
      if i~=1|j~=STEUV(9) error(FT14); end;
      if i1~=1|j1~=STEUV(10) error(FT15); end
      zqw(2:end+1)=zqw(1:end); zqw(1)=1; nqw(2:end+1)=nqw(1:end); nqw(1)=1;
   else zqw=1; nqw=1; end
   [i,j]=size(DGG);
   if i~=1|j<3|STEUV(7)~=1 error(FT09)
   else Vz1=DGG(1); k1=DGG(2); Ns1=DGG(3); zs1=DGG(4:end); end
   if length(zs1)~=Ns1|k1<0|rem(k1,1)~=0|Ns1<0|rem(Ns1,1)~=0 error(FT09); end
   if isempty(zs1) zs1=1; end
elseif nargin==5 [afu,ast]=deal(varargin{:}); Vz1=1;k1=0;Ns1=0;zs1=1;zqw=1;nqw=1; 
elseif nargin==4 afu=deal(varargin{1}); ast=1;Vz1=1;k1=0;Ns1=0;zs1=1;zqw=1;nqw=1; 
else afu=1; ast=1; Vz1=1; k1=0; Ns1=0; zs1=1; zqw=1; nqw=1;
end;
if MUEF==0
   if isempty(afu) afu=1; end; if isempty(ast) ast=1; end
   [i,j]=size(ast); if i~=1|(STEUV(5)<0&(j~=-STEUV(5))) error(FT10); end
   [i,j]=size(afu); if i~=1|(STEUV(4)<0&(j~=-STEUV(4))) error(FT11); end
   if afu(1)~=1 afu(2:end+1)=afu(1:end); afu(1)=1; end
   if ast(1)~=1 ast(2:end+1)=ast(1:end); ast(1)=1; end; end
[znum,snum]=size(n); [zden,sden]=size(ns2); eps1=1e-14; % da eps zu klein ist!
if ((znum~=1)|(zden~=1)) error(FT12); end
if ((snum>9)|(sden>9)) disp(WT01); end
if abs(ns2(sden))<eps1 i=1; j=[1 0];
   while abs(ns2(sden-i))<eps1 i=i+1; j(i+1)=0; end; ns2=deconv(ns2,j);
   sden=sden-i; end
if abs(n(snum))<eps1 i=1; j=[1 0];
   while abs(n(snum-i))<eps1 i=i+1; j(i+1)=0; end; n=deconv(n,j);
   snum=snum-i; end
if ns2(1)~=1 error(FT13); end;
j=0; while abs(n(j+1))<eps1; j=j+1; end; Vz2=n(j+1);
     for i=(1+j):snum   n(i)=n(i)/Vz2;   end; k2=j;
if k2==0 error(FT03); end; if (k2>3) disp(WT02); end
%****************** Vorbereitung der Optimierung *****************************
zs2=1; for i=(2+j):snum zs2(i-j)=n(i); end; Ns2=snum-1-j; M=sden-1;
K=[1 -1 0 0 0 0;1 -2 1 0 0 0;1 -3 3 -1 0 0;1 -4 6 -4 1 0;1 -5 10 -10 5 -1];
if MUEF~=0 fos=input('Optimization: undelayed(any)/delayed(1) =');
else fos=STEUV(1); end
  if isempty(fos)|fos~=1 fos=0;
  else if STEUV(6)~=2&STEUV(6)~=1 disp(KT01); end; end;
if STEUV(6)~=2&STEUV(6)~=1 disp(KT02); end
T='AT00'; AT01='Number integral parts';AT02='Signal wave form';AT03='Polynomial';
if MUEF~=0 ans=input([AT01 '(Enter or 1), ' AT02 '(2), ' AT03 '(3): ']);
else ans=STEUV(8);
     if ans>1 ans=ans+1; elseif ans<1 ans=ans+2; end; %Anpassung: STEUV/Menü
     if STEUV(6)~=2&STEUV(6)~=1 T(4)=48+ans; disp(eval(setstr(T))); end
end
if isempty(ans) ans=1; end
if ans==1 Zqw=0;
   if MUEF~=0 I=input('Order of the controller integral part (Enter=1) = ');
   else I=STEUV(9); end; if isempty(I) I=1; end
   if ((I<0)|(I>5)|(rem(I,1)~=0)) error(FT02); end  %keine sinnvolle Eingabe
   Nw=0; nw=1; Nqw=I; for i=2:(I+1) nqw(i)=K(I,i); end
elseif ans==2
   if MUEF~=0 I=input('Step(1), Ramp(2), Parabola(3): ');
   else I=STEUV(9)+1; end %"+1": Anpassung STEUV an Menü
   if     I==1 Zqw=0;                %zqw=1 durch Initialisierung schon gesetzt
   elseif I==2 Zqw=0;                %statt Zqw=1 und zqw=T*[0 1]; >>>deconv!<<<
          if STEUV(6)~=2 disp(WT09); end %Verstärkungsfaktoren + Sollwert-  
   elseif I==3 Zqw=1; zqw=[1 1];         %Totzeit gehen nicht in Entwurf ein
          if STEUV(6)~=2 disp(WT10); end %statt Zqw=2 und zqw=T^2*[0 0.5 0.5];
   else   if MUEF==1 error(FT01); else error(FT18); end; end
   Nw=0; nw=1; Nqw=I; for i=2:(I+1) nqw(i)=K(I,i); end
elseif ans==3
   if MUEF~=0 Zqw=input('Numerator order= '); T='Input of z0= ';
      for i=1:Zqw T(14)=48+i; T=setstr(T); zqw(1+i)=input(T); end
              Nqw=input('Denominator order= '); T='Input of n0= ';
      for i=1:Nqw T(14)=48+i; T=setstr(T); nqw(1+i)=input(T); end
   else Zqw=STEUV(9); Nqw=STEUV(10);%zqw durch varargin4 schon verfügbar
   end                              %nqw durch varargin5 schon verfügbar
   I=0; j=0; hiv=nqw; nw=nqw;
   while j<eps [hiv,r]=deconv(hiv,[1 -1]); j=sum(r);  % "~j" oder "j==0" wegen
      if j<eps I=I+1; nw=hiv; end; end; Nw=Nqw-I;     %Ungenauigkeiten unzulässig
else if MUEF==1 error(FT01); else error(FT17); end; end
pz=roots(ns2); Is=0;
for i=1:M
    if pz(i)==1 Is=Is+1; end; end
if Is>0 if STEUV(6)~=2&STEUV(6)~=1 disp(KT03); disp(Is); end
   if I>0
    if MUEF~=0 ans=input('Application for controller design? no(0)/yes(number)');
      else ans=STEUV(11); end
      if ans~=0
         if ((rem(ans,1)~=0)|(ans>Is)|(ans<0))
            if MUEF==1 error(FT02); else error(FT19); end; end; hiv=I; I=I-ans;
         if I<0 I=0; if STEUV(6)~=2&STEUV(6)~=1 disp(KT04); disp(hiv); end; end
         if fos==1&STEUV(6)~=2 disp(WT03); end
end; end; end
DT01='"DEE" between controller and '; DT02=' insert? no(any)/yes(1): ';
if MUEF~=0 if fos==0 ans=input([DT01 'continuous plant' DT02]);
           else      ans=input([DT01 'disturbance variable action' DT02]); end
else ans=STEUV(7); end; if isempty(ans)|ans~=1 ans=0; end
if ans==1 if STEUV(6)~=2&STEUV(6)~=1 disp(KT05); end
   if MUEF~=0 
     Vz1=input('Input of gain Vz1= ');
     k1 =input('Input of dead time k1= ');
     Ns1=input('Input of polynomial order= ');
     if (isempty(Vz1)|isempty(k1)|isempty(Ns1)|(rem(Ns1,1)~=0)) error(FT01); end
     T='Input of g0= ';
     for i=1:Ns1 T(14)=48+i; T=setstr(T); zs1(1+i)=input(T); end; end
end %Für MUEF=0 stehen alle Variablen durch varargin3 schon zur Verfügung
k=k1+k2; Vz=Vz1*Vz2;
if fos==1 %*********** Beginn Störoptimierung ***********************************
Ns2str=0; zs2str=1;
if MUEF~=0 ans=input(IT01);
else ans=STEUV(3)+1; end %"+1" zur Anpassung von STEUV-Def. an Input-Abfrage
  if isempty(ans)|ans~=1 ans=0; end
  if ans==1 Ns2ste=Ns2; zs2ste=zs2;
  else
    if Ns2==0 error(FT04); end
    nz=roots(zs2); l=Ns2;
    if STEUV(3)==-1 STEUV(3)=l; end %Vollkompensation = max. Anzahl klein-'el'
    nz=sort(nz);          %wesentlich falls nicht vollständig kompensiert wird
    if STEUV(6)~=2&STEUV(6)~=1 disp(KT06); disp(nz); end; Ns2ste=0; zs2ste=1;
    while l~=0
       if imag(nz(l))==0
          if STEUV(6)~=2&STEUV(6)~=1 disp([KT07 ': ']); disp(nz(l)); end; anzn=1;
       else
         if STEUV(6)~=2&STEUV(6)~=1 disp([KT07 'n: ']);disp(nz(l))
                                                       disp(nz(l-1));end; anzn=2;
       end
       if MUEF~=0 ans=input('yes(any)/no(1)= ');
          if isempty(ans)|ans~=1 ans=0; end
       else ans=STEUV(3);
            if ans~=0 ans=ans-anzn; 
               if ans<0 STEUV(3)=0; if STEUV(6)~=2 disp(WT04); end
               else  STEUV(3)=ans;
               end; ans=0; 
            else ans=1; end
       end
       if MUEF==0&STEUV(6)~=2&STEUV(6)~=1 
          if ans==1 disp('No'); else disp('Yes'); end; end
       for i=1:anzn
           if ans~=1
              Ns2str=Ns2str+1; zs2str(Ns2str+1)=0; arbvek=1; for j=1:Ns2str
                               arbvek(j+1)=zs2str(j+1)-zs2str(j)*nz(l+1-i);
                               end; zs2str=arbvek;
           else
              Ns2ste=Ns2ste+1; zs2ste(Ns2ste+1)=0; arbvek=1; for j=1:Ns2ste
                               arbvek(j+1)=zs2ste(j+1)-zs2ste(j)*nz(l+1-i);
                               end; zs2ste=arbvek;  end
       end; l=l-anzn;
     end %"while" ist sinnvoll da komplexe Nullstellen im Paar kompensieren!
  if Ns2ste~=0 if STEUV(6)~=2 disp(WT05); end; end; end
mst=k1+2*k2+Ns1+Ns2+Ns2ste+Zqw;
if STEUV(6)~=2&STEUV(6)~=1 disp(KT08); disp(mst); end
if MUEF~=0
Ast=input(IT02);
     if isempty(Ast) Ast=0; end
else Ast=abs(STEUV(5));
end
 if ((Ast<0)|(Ast>4)|(rem(Ast,1)~=0))
    if MUEF==1 error(FT02); else error(FT20); end; end
 if Ast~=0 
    if MUEF~=0
      ans=input(IT03);
      if isempty(ans)|ans~=1 ans=0; end
    else if STEUV(5)>0 ans=0; else ans=1; end; end
    if ans==1&MUEF==1  %Menüführung und externe Eingabe
       T='Input of a0 =';
       for i=1:Ast T(14)=48+i; T=setstr(T); ast(i+1)=input(T); end
    elseif ans==0      %Maschinelle Berechnung
     yg=ones(1,Ast+1); Er=Ns2ste+Zqw-Nqw;
     if Er>0 yg(Ast+1+Er)=0; end
     a1=deconv(conv(yg,nqw),conv(zs2ste,zqw));
     if STEUV(6)~=2&STEUV(6)~=1 disp(KT09); end
     for i=1:Ast ast(1+i)=a1(1+i); end; % wegen Polynomdivision nötig
    end; %ans==1&MUEF==0: Werte durch varargin2 schon zur Verfügung
 end 
AV=Nqw-1+M+Ast; BV=k+Ns1+Ns2ste-1+Ast;
cp=conv(nqw,ns2);     CP=Nqw+M;      % mit cp0=1, CP echte Ordnung
dp=conv(zs1,zs2ste);  DP=Ns1+Ns2ste; %  "  dp0=1, DP   "      "
if BV>0
   hiv=-1; cpe=-cp;
     if BV>1  hiv(BV)=0; cpe(BV+CP)=0; end
   ma=toeplitz(cpe,hiv); end
hiv=-conv(dp,ast); HIV=DP+Ast;        % mit hiv0=1
for i=0:HIV
    ma(k+i,BV+1)=hiv(i+1); end
for j=0:(AV-Ast-1)
    for i=0:DP
        ma(k+i+j+1+Ast,BV+2+j)=-dp(i+1); end; end
hiv=0; hiv(BV+CP)=0;
for i=1:CP
    hiv(i)=cp(i+1); end
erg=inv(ma)*hiv',
bv=1; if BV>0
         for i=1:BV
             bv(1+i)=erg(i); end; end
av=erg(1+BV);                        % entspricht V0
for i=1:Ast
    av(1+i)=ast(1+i)*av(1); end
for i=1:(AV-Ast)
    av(1+Ast+i)=erg(1+BV+i); end
out1=av/Vz;                          % Reglerzähler mit einmultipliziertem VR
out2=conv(conv(bv,zs2str),nw);       % I-Anteile werden später zugerechnet
                                     % nur für ROSI: ver=BV+Ns2str+Nw;
end
%************ Führungsoptimierung **********************************************
Ms2str=0; ns2str=1;
if MUEF~=0 
  if fos==0 ans=input(IT04);
  else      ans=input(IT05); end
else ans=STEUV(2)+1; end %"+1" zur Anpassung von STEUV-Def. an Input-Abfrage
  if isempty(ans)|ans~=1 ans=0; end
  if ans==1 Ms2ste=M; ns2ste=ns2;
  else
    if M==0 error(FT05); end
    pz=roots(ns2); l=M;
    if STEUV(2)==-1 STEUV(2)=l; end   %Vollkompensation = max. Anzahl klein-'el'
    pz=sort(pz);            %wesentlich falls nicht vollständig kompensiert wird
    if STEUV(6)~=2&STEUV(6)~=1 disp(KT10); disp(pz); end; Ms2ste=0; ns2ste=1;
    while l~=0
       if imag(pz(l))==0
          if STEUV(6)~=2&STEUV(6)~=1 disp([KT11 'pole: ']); disp(pz(l));
          end;    anzn=1;
       else
          if STEUV(6)~=2&STEUV(6)~=1 disp([KT11 'poles: ']); disp(pz(l))
                                                             disp(pz(l-1));
          end;    anzn=2;
       end
       if MUEF~=0 ans=input('yes(any)/no(1)= ');
          if isempty(ans)|ans~=1 ans=0; end
       else ans=STEUV(2);
            if ans~=0 ans=ans-anzn; 
               if ans<0 STEUV(2)=0; if STEUV(6)~=2 disp(WT06); end
               else  STEUV(2)=ans;
               end; ans=0; 
            else ans=1; end   
       end
       if MUEF==0&STEUV(6)~=2&STEUV(6)~=1 
          if ans==1 disp('no'); else disp('yes'); end; end
       for i=1:anzn
           if ans~=1
              Ms2str=Ms2str+1; ns2str(Ms2str+1)=0; arbvek=1; for j=1:Ms2str
                               arbvek(j+1)=ns2str(j+1)-ns2str(j)*pz(l+1-i);
                               end; ns2str=arbvek;
           else
              Ms2ste=Ms2ste+1; ns2ste(Ms2ste+1)=0; arbvek=1; for j=1:Ms2ste
                               arbvek(j+1)=ns2ste(j+1)-ns2ste(j)*pz(l+1-i);
                               end; ns2ste=arbvek;  end
       end; l=l-anzn;   
     end %"while" ist sinnvoll: komplexe Pole im Paar kompensieren!
  if Ms2ste~=0 if STEUV(6)~=2 if fos==0 disp(WT07); else disp(WT08); end;end;end;
end
mfu=k+Ns1+Ns2+Ms2ste+Zqw;
if STEUV(6)~=2&STEUV(6)~=1 disp(KT12); disp(mfu); end
if MUEF~=0
Afu=input(IT06);
 if isempty(Afu) Afu=0; end
else Afu=abs(STEUV(4)); end 
 if ((Afu<0)|(Afu>4)|(rem(Afu,1)~=0)) 
    if MUEF==1 error(FT02); else error(FT21); end; end
 if Afu~=0 
    if MUEF~=0
      ans=input(IT07);
       if isempty(ans)|ans~=1 ans=0; end
    else if STEUV(4)>0 ans=0; else ans=1; end; end
    if ans==1&MUEF==1 %Menüführung und externe Eingabe
       T='Input of a0 =';
       for i=1:Afu T(14)=48+i; T=setstr(T); afu(i+1)=input(T); end
    elseif ans==0     %Maschinelle Berechnung
     yg=ones(1,Afu+1); Er=M+Zqw-Nqw;
     if Er>0 yg(Afu+1+Er)=0; end
     a1=deconv(conv(yg,nqw),conv(ns2,zqw));
     if STEUV(6)~=2&STEUV(6)~=1 disp(KT13); end
     for i=1:Afu afu(1+i)=a1(1+i); end; % wegen Polynomdivision nötig
    end; %ans==1&MUEF==0: Werte durch varargin1 schon zur Verfügung
 end 
AVf=Nqw-1+Ms2ste+Afu; BVf=k+Ns1+Ns2-1+Afu;
cp=conv(nqw,ns2ste);  CP=Nqw+Ms2ste;      % mit cp0=1, CP echte Ordnung
dp=conv(zs1,zs2);     DP=Ns1+Ns2;         %  "  dp0=1, DP   "      "
if BVf>0
   hiv=-1; cpe=-cp;
      if BVf>1  hiv(BVf)=0; cpe(BVf+CP)=0; end
   ma=toeplitz(cpe,hiv); end
hiv=-conv(dp,afu); HIV=DP+Afu;        % mit hiv0=1
for i=0:HIV
    ma(k+i,BVf+1)=hiv(i+1); end
for j=0:(AVf-Afu-1)
    for i=0:DP
        ma(k+i+j+1+Afu,BVf+2+j)=-dp(i+1); end; end
hiv=0; hiv(BVf+CP)=0;
for i=1:CP
    hiv(i)=cp(i+1); end
erg=inv(ma)*hiv';    %FALSE for Afu==0 and hiv=-1
if Afu==0&hiv==-1 erg=1; end %PATCH ONLY: Afu and afu are OK,
                             %patch for case: 1VH-A [Afu Ast], Afu=0, Ast>0
%erg should increase with Afu, but beeing independent on Ast
bvf=1; if BVf>0
         for i=1:BVf
             bvf(1+i)=erg(i); end; end
avf=erg(1+BVf);                      % entspricht V0
for i=1:Afu
    avf(1+i)=afu(1+i)*avf(1); end
for i=1:(AVf-Afu)
    avf(1+Afu+i)=erg(1+BVf+i); end
out1f=conv(ns2str,avf/Vz);           % Reglerzähler mit einmultipliziertem VR
out2f=conv(bvf,nw);                  % I-Anteile werden später zugerechnet
                                     % nur für ROSI: verf=BVf+Nw;
%************** Ausgabe der Ergebnisse ***************************************
if fos==1 % Filterberechnung
     out3=conv(zs2str,avf)/av(1); out4=out1/out1(1);
else out3=1; out4=1; out1=out1f; out2=out2f; % nur für ROSI: ver=verf         
end
out2=conv(out2,nqw); % Ergänzung um Integralanteile
if LTI==1 out1=tf(out1,out2,inp1,'Variable','z^-1'); %Eigenschaften wie Strecken-
          out2=tf(out3,out4,inp1,'Variable','z^-1'); %LTI-Darstellung, z.B.: T
          out3=[]; out4=[];
else      [j,i1]=size(out1); [j,i2]=size(out2);
          if i1>i2 out2(i1)=0; elseif i1<i2 out1(i2)=0; end
end
