function  [out1,out2,out3,out4] = habod(n,d)
%HABOD    Hilfsfile fuer caspa.m
% Analog zu bod: Berechnung eines betragsoptimalen Reglers ausgehend von
% einer Regelstrecke im Z-Bereich
%	 [out1,out2(,out3,out4)] = bod(num,den)
%                            Fuehrende Nullen in "num" und nachgestellte
%                            Nullen in "den" werden ausgewertet! den(1)=1,
% 	  out1 - Reglerzaehlerpolynom (fallende Z-Potenzen)
%         out2 - Reglernennerpolynom  (    "         "    )
%         out3 - Filterzaehler
%         out4 - Filternennerpolynom  (    "         "    )
%   ==>  (,out3,out4)  nur sinnvoll bei Stoeroptimierung

% SCHLUESSELWOERTER: Hilfsfile

% Vorraussetzung ist das Vorhandensein der m-Files havobo.m

% Copyright (c) 1994 Hagen Jahn
%                    TU Dresden; Elektrotechnisches Institut
%                    Mommsenstr. 13; D-01062 Dresden; Germany

if ((nargin<2)|(nargin==3)|(nargin>4)) error('Eingangsargumentfehler'); end
[znum,snum]=size(n); [zden,sden]=size(d);
eps1=1e-14; % da eps zu klein ist !
if ((znum~=1)|(zden~=1))
   error('Streckenkoeffizienten als Vektoren uebergeben');  end;
if ((snum>9)|(sden>9))
   disp('Ist die Ordnung der Eingangsvektoren tatsaechlich so gross ?'); 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('Nenner(1)=1 erwartet!'); end; j=0;
if nargin==2
   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;
else
   if n(1)~=1 error('Zaehler(1)=1 erwartet!'); end; end
if k==0 disp('Warnung: Algorithmus wird auf sprungfaehige Strecke angewendet')
   end
if k>=3 disp('Warnung: Ist die Totzeit tatsaechlich so gross ?'); end
for i=(2+j):snum  siz(1,1+i-j)=n(i); end
siz(1,1)=Vz; siz(1,2)=snum-1-j;
input('Optimierung: Fuehrung(bel.)/Stoerung(1) ='); if ans==1; else ans=0; end
  opa(1)=ans;
input('Polkompensation: nein(bel.)/ja(1) =');     if ans==1; else ans=0; end
  opa(2)=ans; DKo=0; RZKo=1;
  if ans==1
    if opa(1)==1
       disp('Warnung: Keine optimale Loesung moeglich !');  end
    if sden==1 error('Die eingegebene Strecke besitzt keine Pole'); end
    pz=roots(d); l=sden-1;
    disp('Die Strecke besitzt folgende Pole:'); disp(sort(pz))
    hid=1; shid=1; sRZKo=1; % RZKo=1 ist schon gesetzt
    while l~=0
       if imag(pz(l))==0
          disp('Kompensation des Poles: '); disp(pz(l))
          anzp=1;
       else
          disp('Kompensation der Pole: ');  disp(pz(l))
                                            disp(pz(l-1))
          anzp=2;
       end
       input('nein(bel.)/ja(1)= ');
       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" ist sinnvoll da komplexe Pole im Paar
                            %kompensiert werden sollen
  sden=shid; d=hid; DKo=sRZKo-1; end %RZKo liegt berechnet vor
for i=2:sden    siz(2,1+i)=d(i);  end
siz(2,1)=k; siz(2,2)=sden-1;
input('Reglerintegralanteil: ja(bel.)/nein(1) =');if ans==1 ans=0;else ans=1;end
  opa(4)=ans;
  if ((ans==0)&(opa(1)==1))
disp('Warnung: Vorhandensein eines I-Anteiles vor Stoermischstelle ueberpruefen!')
  end
  if ((ans==0)&(siz(2,2)==0)) error('Keine Berechnung moeglich'); end
input('Ordnung des zu berechnenden Reglerzaehleranteils (ET=0) =');
  if isempty(ans) ans=0; end
  D=ans; b1=0;
  if ((ans==0)&(opa(1)==1))
  error('Nicht sinnvoll - zwei Freiheitsgrade zur Stoeroptimierung notwendig')
  end
    if((ans==0)&(opa(1)==0)&(DKo>1))
      input('"gedaempfter" D-Anteil: nein(bel.)/ja(1) =');
        if ans==1   disp('Warnung: Nur in Sonderfaellen brauchbar!!');
        else        ans=0;  end
      opa(5)=ans;
    else opa(5)=0; end
opa(3)=D; RZO=DKo+D;
if RZO>4 disp('Warnung: Ist diese hohe Reglerordnung sinnvoll ?'); end
[R]=havobo(siz,opa);
if opa(1)==0      %Fuehrungsoptimierung
   if opa(5)==0   %ohne Daempfung des D-Anteiles
      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           %Mit Daempfung des D-Anteiles
                  %Es war D=0 vorgegeben worden!
   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);
   pz=roots([1 (K2+K0)/K3 1+(K1/K3) K0/K3]);
   disp('Die kubische Gleichung fuer b1 liefert folgende Loesungen:')
   disp(pz)
   input('Auswahl einer Loesung (1,2,3): ')
   if (ans<1)|(ans>3) error('Falsche Eingabe')
   elseif imag(pz(ans)~=0) error('Komplexe Reglerparameter sind unzulaessig')
   elseif (pz(ans)<=-1)|(pz(ans)>0)
         disp('Warnung: nur  -1 < b1 < 0  entspricht Daempfungszeitkonstante')
   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   % Stoe”roptimierung
   if D==1
      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('Keine Berechnung moeglich - zu niedrige Streckenstruktur');end
      pe=pzae/pnen; qu=(R(1,1)*R(2,3)-R(2,1)*R(1,3))/pnen;
      V0=roots([1 pe qu])';
disp('Loesungen fuer die Verstaerkung des offenen Kreises: (V01 V02)');disp(V0)
      input('Auswahl: V02(bel.)/V01(1)= ');
      if ans==1  erg=V0(1); else  erg=V0(2);  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
      disp('Der Loesungsalgorithmus benoetigt Anfangswerte !')
      input('Vorgabe: maschinell(bel.)/von Hand(1) =');
      if ans==1
      input('Vorgabe einer: PID-Reglereinstellung(bel.)/Reglernullstelle(1) =');
              if ans==1
                 input('Reglernullstelle =');   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
                 input('Anfangswert fuer VR= '); anfv(1)=ans*siz(1,1);
                 input('Anfangswert fuer d1= '); anfv(2)=ans*anfv(1);
                 input('Anfangswert fuer d2= '); anfv(3)=ans*anfv(1); j=0; end
      else    spz=size(d)-1;
         if spz==0  error('Die Strecke besitzt keine Pole (mehr)'); end
         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]=havobo(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]);
disp('Anfangswert fuer Verstaerkung des offenen Kreises: (V01 V02)'); disp(V0)
      input('Auswahl: V02(bel.)/V01(1) =');
      if ans==1  erg=V0(1); else  erg=V0(2);  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);
      disp('Anfangswerte des Suchalgorithmus (VR d1 d2) :')
      disp([anfv(1)/siz(1,1) anfv(2)/anfv(1) anfv(3)/anfv(1)])
      end;  erg=fsolve('bodst3',anfv',[],[],R); erg=erg'/siz(1,1);
			% Nutzng von "fsolve" vom 21.11.91 2496byte
   else  % "D" wurde nicht mit "1" oder "2" vorgegeben
      error('Dieser Reglergrad ist nicht vorbereitet')
   end
   rzs=sum(erg);
end
out1=conv(erg,RZKo);
if opa(1)==1  % Filterberechnung
   [j,i1]=size(erg); out3(i1)=0;
   out3(1)=rzs/out1(1);
   out4   =erg/out1(1);
else out3=1; out4=1; end
if    ((-opa(4)*b1)~=0)  out2=[1 b1-opa(4) -opa(4)*b1];
elseif ((b1-opa(4))~=0)  out2=[1 b1-opa(4)];
else                     out2=[1];  end
[j,i1]=size(out1); [j,i2]=size(out2);
if i1>i2 out2(i1)=0; elseif i1<i2 out1(i2)=0; end
end
