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

% KEY WORD: supporting file

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

if ~LL                           %case: PID-family controller
   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 ***
   rz=x(1:(D+1))';               %D order controller numerator;  D=0: gain only
   rn=x((D+1):(TAZ))';           %B order controller denominator;B=0: minimum form
   % ***
   rn(1)=1;                      %b0=1; Integral part extra!
   if IP rn=conv(rn,[1 -1]); end %with integral part
else                             %case: Lead/Lag controller
   TAZ=D+B;                      %no explicit gain
   if D==2 TAZ=TAZ-1; end        %second order only with dof=3
   q=zeros(TAZ,1);               %preparation of the fsolve q(ii) computation
   alpha=x(1); beta=x(2);        %definition according paper/conclusion by analogy
   if D==1  rz=[alpha  -(alpha-1)]; %first order Zanasi: identically with
            rn=[beta   -(beta-1) ]; %first order approximation of differentiation
   else     gamma=x(3);          %second order: three unknown quantities (dof=3)
            rz=[alpha  -(alpha-1+gamma) gamma]; %first+second order diff.
            rn=[beta   -(beta-1+gamma+(beta-alpha)/3)  (gamma+(beta-alpha)/3)];
   end
end

if ~PF
   [zz,nz,T]=tfdata(SZ,'v');  %concerning sampling time T
   CZ=tf(rz,rn,T,'variable','z^-1'); %packed controller in z-domain 
   OL=series(CZ,SZ);          %open loop
   SysFB=tf(1,1,T);           %feedback transfer function
   CL=feedback(OL,SysFB);     %closed loop
   [a,b,Tcl]=tfdata(CL,'v');  %polynomials: closed loop; 'v': otherwise cell array
else                          % * no Lead/Lag already checked in calling function
   zfi=sum(rz(1:D+1))/rz(1);  %Filter numerator according final value theorem
   Vr=rz(1);                  %controller gain
   [zz,nz,T]=tfdata(SZ,'v');%plant-polynomials:numerator/denominator;sample time
   a=conv(conv(zz,zfi),Vr); %Compensat.: ControllerNumerator=FilterDenominator
        CZ=tf(rz,rn,T,'variable','z^-1'); %see above - packed controller in z-domain 
        OL=series(CZ,SZ);          %see above - open loop
        SysFB=tf(1,1,T);           %see above - feedback transfer function
        CL=feedback(OL,SysFB);     %see above - closed loop
        [CLnumOF,b,Tcl]=tfdata(CL,'v');  %see above/closed loop without filter
                                   %same denominator with/ without filter
end % a and b are available (numerator / denominator polynomials closed loop)

%check of computability
  oz=length(a); on=length(b);   % order of numerator and denominator
  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 Magnitude 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
