J=bod_esp(p,a_o,b_o,c_o,d_o,EZW,QQ,BMEA,ZVe,NVe,NV,TOL)

function J=bod_esp(p,a_o,b_o,c_o,d_o,EZW,QQ,BMEA,ZVe,NVe,NV,TOL)
%BOD_ESP  auxiliary file for bod_es.m - optimize structure in z-domain
%                J = bod_esp(p,a_o,b_o,c_o,d_o,EZW,QQ,BMEA,ZVe,NVe,NV,TOL)
%         obsolete call:
%         [J,nfun] = bod_esp(p,a_o,b_o,c_o,d_o,EZW,QQ,BMEA,ZVe,NVe,NV)
%         for definition of plant and control structure see file:
%         strucxxx.m (example: struc000.m)
%	  Here: computation of partial derivatives applicable for products of
%		    two unknown only (otherwise see MATLAB manual)

%  KEY WORDS: Digital Amount Optimum, optimization in z-domain,
%             auxiliary file for bod_es

%  Copyright (c) 2011 Dr. Gert-Helge Geitner; TU Dresden, Fak. ET,
%                     Elektrotechnisches Institut (ETI); Mommsenstr. 13;
%                     D-01062 Dresden, Germany;
%                     Tel./Fax ++49-351-463-32917/33655
%                     [e-mail]

%unpack of transfered packed fsolve arguments
 BRF=EZW(1,1); MAU=EZW(1,2); TAZ=EZW(1,3); iS=EZW(1,4);
 knr=EZW(1,5);                       % number of already known results
 for i=1:knr KZV_o(i)=EZW(2,i); end; % KZV_o equal to erg
 [i,j]=size(BMEA);j=j/2;BME(:,1:j)=BMEA(:,1:j);BMA(:,1:j)=BMEA(:,j+1:2*j);

KZV_o(MAU)=0;	       % "remaining" coefficients firstly fill up by zeros
for i=1:TAZ
    KZV_o(knr+i)=p(i); % fsolve intermediate results; KZV_o: original value
end
J(TAZ,TAZ)=0;          % set matrix order for results
KZV_s=[1 -1 2 -2 3 -3];% specification
for u=1:TAZ            % TAZ = number of partial derivatives per equation
    q=zeros(6,TAZ);    % preparation of fsolve q(ii) computation
% manual pp.3-87 wrong Jacobian matrix design
    KZV=KZV_o;
    for w=1:6          % f(x=1);f(x=-1);f(x=2);f(x=-2);f(x=3);f(x=-3);
        KZV(knr+u)=KZV_s(w);
 	a=a_o; b=b_o; c=c_o; d=d_o;
% START OF INNER LOOP-ALGORITHM
% unpack name vector of controler parameters and assign numerical values:
 i=0; n0=0; k=length(NV);
 while i<k n0=n0+1; i=i+1; j=1; vtext=0; %increase "i" by "2" because of: "; "
       while NV(i)~=';' vtext(j)=NV(i); i=i+1; j=j+1; end
       eval([vtext '=KZV(' num2str(n0) ');']); i=i+1; end

% unpack numerator vector of particular blocks each:
 i=0; n1=0; k=length(ZVe);
 while i<k n1=n1+1; i=i+1; j=1; vtext=0;
       while ZVe(i)~=';'  vtext(j)=ZVe(i); i=i+1;  j=j+1; end
       eval(['rz' num2str(n1) '=' vtext ';']); end

% unpack denominator vector of particular blocks each:
 i=0; n=0; k=length(NVe);
 while i<k n=n+1; i=i+1; j=1; vtext=0;
       while NVe(i)~=';'  vtext(j)=NVe(i); i=i+1;  j=j+1; end
       eval(['rn' num2str(n) '=' vtext ';']); end

% generate unconnected diagonal system:
for i=1:n
    RZ=eval(['rz' num2str(i)]);RN=eval(['rn' num2str(i)]);
    zl=length(RZ); nl=length(RN); if nl<zl RN(zl)=0; end    % because of tf2ss
    [at,bt,ct,dt]=tf2ss(RZ,RN);
    [a,b,c,d]=append(a,b,c,d,at,bt,ct,dt); end

% compute connected system:
 [ac,bc,cc,dc]=connect(a,b,c,d,QQ,BME(iS,BRF),BMA(iS,BRF));
 [acm,bcm,ccm,dcm]=minreal(ac,bc,cc,dc,TOL);  %===> unconditional necessary!!!
 [a,b]=ss2tf(acm,bcm,ccm,dcm);

% check of computability:
  oz=length(a); on=length(b);   % order of numerator and denominator
  o=max(oz,on);     % protection against denominator polynomial approximations
  if oz>on b(oz)=0; elseif on>oz a(on)=0; end
  Z=TAZ;			% number of rows for Matrix K
  S=o-1;			% number of columns for Matrix K
  if TAZ>=o disp('Warning: order < number of unknown!')
     a(TAZ+1)=0; b(TAZ+1)=0; o=TAZ+1;
  end

% if u==1&w==1        % S(f(x=1)); old fragment: compute matrix K once only
for i=1:Z                       %Computation of matrix K: columns=S; rows=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; % end

for ii=1:TAZ                                        %start of BOD algorithm 
    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(w,ii) = q(w,ii) + x * K(ii,i);  end; end  % end of BOD algorithm
%END OF INNER LOOP-ALGORITHM
end  %end of for-loop: "w=1:6" to use KZV
%Jacobian computation:
D(1,:)=q(1,:)-q(2,:); D(2,:)=q(1,:)+q(2,:);
D(3,:)=q(3,:)-q(4,:); D(4,:)=q(3,:)+q(4,:); D(6,:)=q(5,:)+q(6,:);
D(5,:)=D(6,:)-D(2,:); D(6,:)=D(6,:)-D(4,:);
PK(1,:)=(8*D(1,:)-D(3,:))/12; PK(2,:)=13/24*D(5,:)-2/3*D(6,:);
PK(3,:)=(D(3,:)-2*D(1,:))/4;  PK(4,:)=(D(5,:)-1.6*D(6,:))/92;
J(u,:)=PK(1,:)+PK(2,:)*KZV_o(knr+u)+...
       PK(3,:)*KZV_o(knr+u)^2+PK(4,:)*KZV_o(knr+u)^3;
end   %end of for-loop: "TAZ partial derivatives"
%for obsolete fsolve call a second output argument was necessary: nfun
%nfun=TAZ^2;        %because of definition in MATLAB manual pp.3-87
%J(u,:)=J(u,:)+(2*w-3)*q(1,:);   % -f(x=1)+f(x=2) %obsolete fragment