function R=vobo(siz,opa)
%VOBO    supporting function for bod.m; ATTENTION: do NOT call stand-alone
%        R=vobo(siz,opa)
% Reglereinstellung nach Digitalem Betragsoptimum (BOD); Vorbereitung für
% BOD-Optimierung ausgehend von Strecke im Z-Bereich (siz) und ausgewählter
% Optimierungsart (opa) durch Berechnung der Rx-Koeffizienten.
% R=[R2F R2L0 R2L1 ... R2LD R2Q0 ... R2QD R2ZL R2ZP;
%     .                 .    .        .   R4ZL R4ZP;
%     .                 .   RAQ0 ... RAQD  0    0  ;
%    RAF  .    .   ... RALD  0   ...  0    0    0  ;]
%
% Es werden nur die entsprechend "opa" tatsächlich benötigten Zeilen und
% Spalten von "r" berechnet.
%
% siz = [Vz N n1 ... nN;         opa = [fos kon D I gon]
%        k  M m1 ... mM]
% M  : steht bei Kompensation für M-quer
% fos: führungs- (0) oder störoptimal (1)
% kon: Kompensation (1) oder nicht (0)
% D  : Ordnung des mit "R" zu berechnenden Reglerzähler(anteil)s
% I  : Reglerintegralanteil (1) oder nicht (0)
% gon: gedämpfter D-Anteil (1) oder nicht (0) - nur für Sonderfall D=0
%      und Führungsoptimierung (z.Z.)

% SCHLÜSSELWORTE: Hilfsfile

%  Copyright (c) 2000-2009 Dr. Gert-Helge Geitner; TU Dresden, Fak. ET,
%                     Elektrotechnisches Institut (ETI); Mommsenstr. 13;
%                     D-01062 Dresden, Germany;
%  Literatur zum Digitalen Betragsoptimum: Geitner, G.-H.:
%                     Entwurf digitaler Regler für elektrische Antriebe.
%                     VDE-Verlag GmbH Berlin und Offenbach 1996
%                     ISBN 3-8007-1847-2

% Es werden folgende allgemeine Parameter definiert:
alf=-opa(4); D=opa(3); N=siz(1,2); k=siz(2,1); M=siz(2,2);
I1=max([(N+k)   (M+1-k)]);
I2=max([(N+k-1) (M+2-k)]);
Ks=max([M+2 I1+D N+D I2]);   %Ks={M+1; I+D; M+2; N+D; I2 }
Kz=max([D+1 2]);             %Kz={D;   D+1; 2;    %    % }
K=kmatrix(Kz,Ks);            %Grösste Koeffizienten der K-Matrix

% Die Werte RAF werden für jede Variante benötigt:
maxi=M+1; mini=-1; zei=2;    %RAF: Kz=D; Ks=maxi;
m=vaib(maxi,mini,zei,siz);
for x=1:(D+1+opa(5))
    R(x,1)=0;
    for i=x:(M+1)
        var=0;
           for j=(0+1-mini):((M-i+1)+1-mini)
               var=(alf*m(j-1)+m(j))*(alf*m(j+i-1)+m(j+i))+var;
           end
        R(x,1)=K(x,i)*var+R(x,1);
    end
end
% Die Werte RAL werden für jede Variante benötigt."D" und "opa(5)"
% entscheiden über die Anzahl.
siz(1,2)=1; I=I1;
maxi=N+k+I+(2*D); mini=k-I-1; if mini>0 mini=0; end %RAL: Kz=D+1; Ks=I1+D; 
m=vaib(maxi,mini,zei,siz);
for l=0:D
    for x=1:(D+1+opa(5))
        R(x,l+2)=0;
        for i=x:(I+l)
            var=0;
            for j=(k+1-mini):((N+k)+1-mini)
            o=j-1+mini;
  var=siz(1,2+o-k)*(m(j-i+l)+(alf*m(j-i+l-1))+m(j+i+l)+(alf*m(j+i+l-1)))+var;
            end
            R(x,l+2)=K(x,i)*var+R(x,l+2);
        end
    end
end
siz(1,2)=N;
if opa(1)==1
% Die Werte RAQl werden nur für Störoptimierung benötigt
maxi=N+D; mini=-D; zei=1;    %RAQ1: Ks=maxi;
n=vaib(maxi,mini,zei,siz);
for l=0:D
    for x=1:(D+1)
        R(x,3+D+l)=0;
        for i=x:(N+l)
            var=0;
            for j=(0+1-mini):((N-i+l)+1-mini)
                var=(n(j)*n(j+i-l))+(n(j-l)*n(j+i))+var;
            end
            R(x,3+D+l)=K(x,i)*var+R(x,3+D+l);
        end
    end
end
elseif opa(5)==1
% Die Werte RAZL und RAZP werden nur für verzögerten D-Anteil benötigt.
maxi=2*M+4; mini=-(M+2); zei=2;   %RAZL/RAZP: Kz=2; Ks=maxi/2;
m=vaib(maxi,mini,zei,siz);
for x=1:2
    R(x,3+D)=0;
    for i=x:(M+2)
        var=0;
        for j=(1+1-mini):((M+2)+1-mini)
   var=(alf*m(j-2)+m(j-1))*(m(j-i)+(alf*m(j-i-1))+m(j+i)+(alf*m(j+i-1)))+var;
        end
        R(x,3+D)=K(x,i)*var+R(x,3+D);
    end
end
I=I2; siz(1,2)=1;
maxi=N+k+I-1; mini=k-I-2; if mini>0 mini=0; end  %R-Matrix: Ks=I2;
m=vaib(maxi,mini,zei,siz);
for x=1:2
    R(x,4+D)=0;
    for i=x:I
        var=0;
        for j=(k+1-mini):((N+k)+1-mini)
        o=j-1+mini;
      var=siz(1,2+o-k)*(m(j-i-1)+(alf*m(j-i-2))+m(j+i-1)+(alf*m(j+i-2)))+var;
        end
        R(x,4+D)=K(x,i)*var+R(x,4+D);
    end
end
siz(1,2)=N;
end


function hiv=vaib(maxi,mini,zei,siz)
%Regleroptimierung nach BOD: Vorbereitung für allg. Indexberechnung
%
if (mini<0) hiv(-mini)=0; end
hiv (1-mini)=1;
for i=1:siz(zei,2) hiv(1-mini+i)=siz(zei,2+i); end
hiv(maxi+1-mini)=0; % Vorsicht: Wenn Vektorordnung nicht ohnehin über-
                    % schritten wird, extern maxi:=maxi+1 nötig !!!