function q=bod_rst_eq(x,SZ,D,B,IP)
%bod_rst_eq supporting file of bod_rst.m - structure in z-domain
%        q=bod_rst_eq(x,SZ,D,B,IP)
% **** ATTENTION: calling by bod_rst only!

% KEY WORD: supporting file

%  Copyright (c) 2018. All Rights Reserved.
%                     Dr. Gert-Helge Geitner; TU Dresden, Fak. ET,
%                     Elektrotechnisches Institut (ETI); Mommsenstr. 13;
%                     D-01062 Dresden, Germany;
%                     https://etiema.et.tu-dresden.de/ae2_files/ae_1_0e.htm
TAZ=D+1+B;                    %definition adopted from: zust_au...
q=zeros(TAZ,1);               %preparation of the fsolve q(ii) computation
% *** Attention: x return: column vector ***
R=x(1:(D+1))';                %D order controller R-part; D=0: gain only
R(1)=sqrt(R(1));  %fuer REELLE Loesung muss R(1) POSITIV werden, damit VR
                  %Quadrat gefunden und Ergebnis ueber WRZL bestimmbar
S=x((D+1):(TAZ))';            %B order controller S-part; B=0: minimum form
% ***
S(1)=1;                       %b0=1; Integral part extra!
if IP S=conv(S,[1 -1]); end   % controller S-part with integrator

[zz,nz,T]=tfdata(SZ,'v');      %concerning sampling time T
Rtf=tf(R,1,T,'variable','z^-1'); %packed controller R-part in z-domain 
Stf=tf(1,S,T,'variable','z^-1'); %packed controller S-part in z-domain
OL=series(Stf,SZ);             %open loop
CL=feedback(OL,Rtf);           %closed loop
[a,b,Tcl]=tfdata(CL,'v');      %closed loop polynomials; 'v': or cell array
% NOW a and b available (closed loop numerator / denominator polynomials)

%check of computability
  oz=length(a); on=length(b);  % order of numerator and denominator of CL
  o=max(oz,on);	 % protect in case of denominator polynomial approximations
  if oz>on b(oz)=0; elseif on>oz a(on)=0; end
  Z=TAZ;			% K Matrix: number of rows
  S=o-1;			% K Matrix: number of columns
  if TAZ>=o disp('Warning: order<number_of_unknown_quantities!')
     a(TAZ+1)=0; b(TAZ+1)=0; o=TAZ+1;
  end
%K Matrix calculation of columns=S and rows=Z
for i=1:Z
 for j=1:S
  if     i==1         K(i,j)=j^2;
  else   if j==1      K(i,j)=0;
	     elseif  j==2 K(i,j)=K(i,j-1)+K(i-1,j-1);
	     elseif  j==3 K(i,j)=K(i,j-1)+K(i,j-2)+K(i-1,j-1)+K(i-1,j-2);
	     else      K(i,j)=K(i,j-1)+K(i,j-2)-K(i,j-3)+K(i-1,j-1)+K(i-1,j-2);
         end; 
  end; 
 end;
end
%Discontinuous Amplitude Optimum: general optimization equation system
for ii=1:TAZ
    for i=ii:(o-1), x=0;
        for j=1:(o-i)
            x = x + a(j) * a(j+i) - b(j) * b(j+i);
        end
        q(ii) = q(ii) + x * K(ii,i);
    end;
end
