function [out1,out2,out3,out4] = bod(inp1,inp2,inp3)
%BOD     Calculation of a discontinuous controller using the
%        Discontinuous Magnitude Optimum (BOD) - starting via a plant
%        given in z-domain, typically applying function trans of the BOD toolbox.
%        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)] = bod(inp1,inp2(,inp3))
%                       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)] = bod(inp1(,inp2))
%                       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(1) - any/1 = Optimization regarding: undelayed / delayed inputs
%       inpx(2) - any/1 = pole compensation: no/yes
%       inpx(3) - any/1 = controller integral part: yes/no
%       inpx(4) = number of poles to be compensated (complex: pole pairs;
%                                                  real: poles - and sequence via 
%                 absolute value: descending undelayed / ascending delayed inputs
%       inpx(5) = order of the controller numerator part to be >>calculated<<
%                 (total numerator order: compensation + calculation)
%       inpx(6) - any/1/2 = elimination of text outputs: no / yes(C) / yes(W+C)
%                 (W - warnings; C - comments)
%       inpx(7) - any/1 = implicate special case damped derivative part: no/yes
%       inpx(8) - 0/1/2 - choice from 3 cubic equa. solutions as to b1; inpx(7)=1
%       inpx(9) - any/1 = solution via quadratic equation - choice: V02/V01
%                 (rule of thumb: V02 as the gain of the open loop)
%       inpx(10)- any/1/2 = solution via search algorithm - PID initial values:
%                 automatic / manual control parameters / manual control zeros
%       inpx(11)...inpx(13) = without meaning if inpx(10)=any
%                           = [VR d1 d2]      if inpx(10)=1
%                           = [zero % %]      if inpx(10)=2
%   any == any other input, e.g. enter or 0, matches preferential input setting;
%   >>>Text output eliminations are recommended after getting familiar only!<<<
%   Defaults and completions of missing elements of inpx: zero

% KEY WORDS: Discontinuous Magnitude Optimum, discontinuous controller of
%            PID-type family, control loop, optimization regarding undelayed
%		     input signals, optimization regarding delayed input signals,
%            pre-filter design, optionally pole compensation,
%		     optionally menu driven, optionally LTI syntax

%  ***** Preconditions: existence of bodg, vobo - see BOD toolbox ****
%  Copyright (c) 1997-2023. All Rights Reserved. 
%                     Dr. Gert-Helge Geitner, TU Dresden,
%                     https://wwwpub.zih.tu-dresden.de/~geitner/ae_1_0e.htm
%  Reference regarding the application of the Discontinuous Magnitude Optimum:
%                     Geitner, G.-H.:
%                     "Entwurf digitaler Regler für elektrische Antriebe",
%                     publisher: VDE-Verlag GmbH Berlin und Offenbach 1996
%                     ISBN 3-8007-1847-2

FT01='Input argument error';
FT02='Old syntax: input of the coefficients of the plant as vectors';
FT03='Denominator(1)=1 expected!';
FT04='Input argument control_vector(inpx) - invalid dimension';
FT05='Input argument control_vector(inpx) - element 4 or/and 5 incorrect!';
FT06='Pole compensation impossible: that plant features no poles';
FT07='Calculation impossible with that setting!';
FT08='Not helpful - two degrees of freedom necessary for delayed input signals';
FT09='Incorrect input for a choice of a solution regarding b1';
FT10='Choice of a solution for b1: complex controller parameters are incorrect';
FT11='Optimization for delayed inputs impossible - too less plant structure';
FT12='Delayed input signal (PID): not enough poles for automatic initial values';
FT13='Delayed input signals: controller order > PID not prepared!';
WT01='Warning: order of the input vectors > 9 in fact?';
WT02='Warning: use of the algorithm for systems with feed through!';
WT03='Warning: dead time > 3 in fact?';
WT04='Warning: this setting enables no optimal solution!';
WT05='Warning: number of poles < number of poles intended to be compensated!';
WT06='Warning: automatic increase of the controller order due to complex poles!';
WT07='Warning: make sure integral part before delayed input mixing point!';
WT08='Warning: with damped derivative part - useful in special cases only!!';
WT09='Warning: such high controller order > 4 useful?';
WT10='Warning: only  -1 < b1 < 0  corresponds with a delay time constant';
WT11='Warning: Abortion! K3==0 and N=k=1 would cause zero compensation.';
KT01='The plant features following poles:';
KT02='Compensation of the pole: ';
KT03='Compensation of the poles: ';
KT04='The cubical equation regarding b1 results in the following solutions:';
KT05='Delayed input signals (PI): solutions of the open loops gain: (V01 V02)';
KT06='Delayed input signals (PID): search algorithm needs initial values!';
KT07='Initial value of the gain of the open loop: (V01 V02)';
KT08='Initial values of the search algorithm (VR d1 d2) :';

if (nargin<1)|(nargin>3)|(nargin==2&isa(inp2,'lti'))|...
   (nargin==3&isa(inp3,'lti')) error(FT01); end
if isa(inp1,'lti')          %LTI model syntax used
   if (nargin==3)|~issiso(inp1)|~isdt(inp1) error(FT01); end
   if ~isa(inp1,'tf') inp1=tf(inp1); end
   [n,d]=tfdata(inp1,'v'); LTI=1;
   if nargin==2 MUEF=0; STEUV=inp2; else MUEF=1; STEUV(6)=0; end   
else                        %"old" numerator / denominator syntax used
   if nargin==1 error(FT01); end
   n=inp1; d=inp2; LTI=0;
   if nargin==3 MUEF=0; STEUV=inp3; else MUEF=1; STEUV(6)=0; end
end
[znum,snum]=size(STEUV); if znum>1|snum>13 error(FT04); end
if snum<13 STEUV(13)=0; end %default addition
if STEUV(4)<0|rem(STEUV(4),1)~=0 error(FT05); end
if STEUV(5)<0|rem(STEUV(5),1)~=0 error(FT05); end
          %for remaining STEUV check input check may be used!
[znum,snum]=size(n); [zden,sden]=size(d); eps1=1e-14; %since "eps" too less
if ((znum~=1)|(zden~=1)) error(FT02);  end;
if STEUV(6)~=2 if ((snum>9)|(sden>9)) disp(WT01); end; end
if abs(d(sden))<eps1 i=1; j=[1 0];
   while abs(d(sden-i))<eps1 i=i+1; j(i+1)=0; end; d=deconv(d,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 d(1)~=1 error(FT03); end
j=0; while abs(n(j+1))<eps1 j=j+1; end; Vz=n(j+1);
     for i=(1+j):snum   n(i)=n(i)/Vz; end; k=j;
if STEUV(6)~=2 if k==0 disp(WT02); end; if k>3 disp(WT03); end; end

for i=(2+j):snum  siz(1,1+i-j)=n(i); end
siz(1,1)=Vz; siz(1,2)=snum-1-j;
if   MUEF~=0 ans=input('Optimization: undelayed(any)/delayed(1) ='); 
else ans=STEUV(1); end
  if isempty(ans)|ans~=1 ans=0; end 
  opa(1)=ans;
if   MUEF~=0 ans=input('Pole compensation: no(any)/yes(1) =');
else ans=STEUV(2); end    
  if isempty(ans)|ans~=1 ans=0; end
  opa(2)=ans; DKo=0; RZKo=1;
  if ans==1
     if opa(1)==1&STEUV(6)~=2 disp(WT04); end %delayed input optimization
  if sden==1 error(FT06); end
    pz=roots(d);
    if opa(1)==0 pz=sort(pz); end %STEUV: difference delayed or undelayed inputs
    l=sden-1;
    if STEUV(6)~=2 if l<STEUV(4) disp(WT05); end; end
    if STEUV(6)~=2&STEUV(6)~=1 disp(KT01); disp(pz); end
    hid=1; shid=1; sRZKo=1; % RZKo=1 is already set
    while l~=0 %while loop for pole compensation
       if imag(pz(l))==0
          if STEUV(6)~=2&STEUV(6)~=1 disp(KT02); disp(pz(l)); end; anzp=1;
       else
          if STEUV(6)~=2&STEUV(6)~=1 disp(KT03); disp(pz(l));
                                         disp(pz(l-1)); end; anzp=2;
       end
       if   MUEF~=0 ans=input('no(any)/yes(1)= ');
            if isempty(ans)|ans~=1 ans=0; end
       else ans=STEUV(4);
            if ans~=0 ans=ans-anzp; 
               if ans<0 STEUV(4)=0; if STEUV(6)~=2 disp(WT06); end
               else  STEUV(4)=ans;
               end; ans=1; end
       end
       if MUEF==0&STEUV(6)~=2&STEUV(6)~=1 
          if ans==0 disp('No'); else disp('Yes'); end; end
       for i=1:anzp
           if ans==1
              sRZKo=sRZKo+1; RZKo(sRZKo)=0; arbvek=1; for j=1:(sRZKo-1)
                               arbvek(j+1)=RZKo(j+1)-RZKo(j)*pz(l+1-i);
                               end; RZKo=arbvek;
           else
              shid=shid+1; hid(shid)=0; arbvek=1; for j=1:(shid-1)
                               arbvek(j+1)=hid(j+1)-hid(j)*pz(l+1-i);
                               end; hid=arbvek;  end
       end; l=l-anzp; 
     end %while: is meaningful since pairwise complex pole compensation necessary
  sden=shid; d=hid; DKo=sRZKo-1; end %RZKo is known here
for i=2:sden    siz(2,1+i)=d(i);  end
siz(2,1)=k; siz(2,2)=sden-1;
if   MUEF~=0 ans=input('Controller integral part: yes(any.)/no(1) =');
else ans=STEUV(3); end
  if isempty(ans)|ans~=1 ans=1; else ans=0; end
  opa(4)=ans;
  if ((ans==0)&(opa(1)==1)) if STEUV(6)~=2 disp(WT07); end; end
  if ((ans==0)&(siz(2,2)==0)) error(FT07); end %computation impossible 
if   MUEF~=0 ans=input('Order of the numerator part to be computed (Enter=0) =');
             if isempty(ans) ans=0; end
else ans=STEUV(5); end
  D=ans; b1=0;
  if ((ans==0)&(opa(1)==1)) error(FT08); end %2 degrees of freedom for undelayed
  if((ans==0)&(opa(1)==0)&(DKo>1))           %inputs mecessary at least
    if MUEF~=0 ans=input('"damped" derivative part: no(any)/yes(1) =');
    else ans=STEUV(7); end
    if isempty(ans)|ans~=1 ans=0;
    else if STEUV(6)~=2 disp(WT08); end; end %useful for special cases only
    opa(5)=ans;
  else opa(5)=0; end
opa(3)=D; RZO=DKo+D;
if RZO>4 disp(WT09); end %high controller order meaningful?
[R]=vobo(siz,opa);
if opa(1)==0             %undelaed input optimization
   if opa(5)==0          %without damping of derivative(D)-part
      for i=1:(D+1)
          svek(i,1)=-R(i,1);
          for j=1:(D+1)
              RM(i,j)=R(i,j+1);
               end; end
      erg=(inv(RM)*svek)'/siz(1,1);
   else %with damping of damping(D)-part - RZ part to be computed is known!
   K0=R(1,1)*R(2,2)-R(1,2)*R(2,1);
   K1=R(1,3)*R(2,2)-R(1,2)*R(2,3);
   K2=R(1,3)*R(2,4)-R(1,4)*R(2,3);
   K3=R(1,1)*R(2,4)-R(2,1)*R(1,4);
   if K3==0 disp(WT11), out1=1; out2=1; out3=1; out4=1;
      return, end %N=k=1: would result in zero compensation
   pz=roots([1 (K2+K0)/K3 1+(K1/K3) K0/K3]); pz=sort(pz); %avoid K3=0!
   if STEUV(6)~=2&STEUV(6)~=1 disp(KT04); disp(pz); end
   if MUEF~=0 ans=input('Choice of a solution (1,2,3): ');
   else ans=STEUV(8)+1; end
   if isempty(ans) ans=1; end        %to block accidentally "enter" input via '1'
   if (ans<1)|(ans>3) error(FT09)      %false input
   elseif imag(pz(ans)~=0) error(FT10) %complex controller parameters illegal
   elseif (pz(ans)<=-1)|(pz(ans)>0) disp(WT10); end; b1=pz(ans);
   erg=-(R(2,1)+b1*(R(2,3)+b1*R(2,1)))/((R(2,2)+b1*R(2,4))*siz(1,1));
   end
else   % delayed input optimization
   if D==1 %delayed input PI controller
      pzae=(R(1,4)-R(1,5))*R(2,1)-(R(2,4)-R(2,5))*R(1,1)+R(1,2)*R(2,3)...
                                                           -R(1,3)*R(2,2);
      pnen=(R(1,4)-R(1,5))*R(2,2)-(R(2,4)-R(2,5))*R(1,2);
      if pnen==0 error(FT11); end       %plant order too less
      pe=pzae/pnen; qu=(R(1,1)*R(2,3)-R(2,1)*R(1,3))/pnen;
      V0=roots([1 pe qu])';
      if STEUV(6)~=2&STEUV(6)~=1 disp(KT05); disp(V0); end %solution open loop
      if MUEF~=0 ans=input('Choice: V02(any)/V01(1)= ');
      else ans=STEUV(9); end                       %rule of thumb: V01
      if isempty(ans)|ans~=1 erg=V0(2); else erg=V0(1); end
      erg(2)=(erg(1)*R(2,2)+R(2,1))/(erg(1)*(R(2,4)-R(2,5))-R(2,3));
      erg=erg/siz(1,1);
   elseif D==2 %delayed input PID controller
      if STEUV(6)~=2&STEUV(6)~=1 disp(KT06); end
      if MUEF~=0 ans=input('Choice: automatic(any)/manual(1) =');
      else ans=STEUV(10); if ans==2 ans=1; end; end %summing-up for MUEF=0
      if isempty(ans)|ans~=1 ans=0; end
      if ans==1
if MUEF~=0 
   ans=input('Choice: PID-Setting(any)/ControllerZero(1) =');
else ans=STEUV(10)-1; end
if isempty(ans)|ans~=1 ans=0; end
        if ans==1
           if MUEF~=0 ans=input('Controller zero = '); else ans=STEUV(11); end
           trz=[1 ans]; j=1; n=conv(n,trz); siz(1,2)=snum;
           for i=3:(snum+2)  siz(1,i)=n(i-1);  end
        else
           if MUEF~=0 ans=input('Initial value VR= '); else ans=STEUV(11); end
           anfv(1)=ans*siz(1,1);
           if MUEF~=0 ans=input('Initial value d1= '); else ans=STEUV(12); end
           anfv(2)=ans*anfv(1);
           if MUEF~=0 ans=input('Initial value d2= '); else ans=STEUV(13); end
           anfv(3)=ans*anfv(1); j=0;
        end
      else    spz=size(d)-1;
         if spz==0  error(FT12); end  %plant without poles (no more)!
         pz=roots(d); kpl=sum(abs(imag(pz)));
         if kpl==0  i1=1;
            for i=1:spz
                if pz(i)>=0
                   pzpos(i1)=pz(i); i1=i1+1;
                end; end
            if     i1>1     trz=[1 -min(pzpos)];
            else   trz=[1 -max(pz)];  end;  j=1;
            d=deconv(d,trz); siz(2,2)=sden-1-j; d(sden)=0;
            for i=3:(sden+1) siz(2,i)=d(i-1); end
         else
            anfv=[siz(1,1) siz(1,1) siz(1,1)]; j=0;
         end
      end
      opa(3)=1; [E]=vobo(siz,opa);
      if j==1
        pzae=(E(1,4)-E(1,5))*E(2,1)-(E(2,4)-E(2,5))*E(1,1)+E(1,2)*E(2,3)...
                                                              -E(1,3)*E(2,2);
        pnen=(E(1,4)-E(1,5))*E(2,2)-(E(2,4)-E(2,5))*E(1,2);
        pe=pzae/pnen; qu=(E(1,1)*E(2,3)-E(2,1)*E(1,3))/pnen;
        V0=roots([1 pe qu]);
        if STEUV(6)~=2&STEUV(6)~=1 disp(KT07); disp(V0); end %initial value (IV) 
        if MUEF~=0 ans=input('Choice: V02(any)/V01(1) =');   %'V' for open loop
           else ans=STEUV(9); end
        if isempty(ans)|ans~=1 erg=V0(2); else erg=V0(1); end
        erg(2)=(erg(1)*E(2,2)+E(2,1))/(erg(1)*(E(2,4)-E(2,5))-E(2,3));
        anfv=conv(erg,trz);
        if STEUV(6)~=2&STEUV(6)~=1 disp(KT08); %IV=[VR d1 d2] of search algorithm
           disp([anfv(1)/siz(1,1) anfv(2)/anfv(1) anfv(3)/anfv(1)]); end
      end
      erg=fsolve('bodg',anfv',[],R);  erg=erg'/siz(1,1);
   else error(FT13);               % "D" was not determined by "1" or "2"
   end; rzs=sum(erg);
end
out1=conv(erg,RZKo);                                %controller numerator
if opa(1)==1 [j,i1]=size(erg);                      %delayed input optimization
             out3(i1)=0; out3(1)=rzs/out1(1);       %filter numerator
             out4=erg/out1(1);                      %filter denominator
else out3=1; out4=1; end                            %undelayed input optimization
if    ((-opa(4)*b1)~=0)  out2=[1 b1-opa(4) -opa(4)*b1]; %controller denominator
elseif ((b1-opa(4))~=0)  out2=[1 b1-opa(4)];
else                     out2=[1];  end
if LTI==1 out1=tf(out1,out2,inp1,'Variable','z^-1'); %features like plants LTI
          out2=tf(out3,out4,inp1,'Variable','z^-1'); %representation, e.g.: T
          out3=[]; out4=[];                          %(sampling time)
else      [j,i1]=size(out1); [j,i2]=size(out2);
          if i1>i2 out2(i1)=0; elseif i1<i2 out1(i2)=0; end
end
