function [nuco,deco,nufil,defil] = casbod(n,d,code)
%CASBOD   Hilfsfunktion fuer casc.m
%         [nuco,deco,nufil,defil] = casbod(n,d,code)
% Berechnet einen betragsoptimalen Regler; Ausgangspunkt: Strecke im z-Bereich

% SCHLUESSELWOERTER: Hilfsfile

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

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

[rnu,cnu]=size(n); [rde,cde]=size(d);
eps1=1e-14; 	% da eps zu klein ist !
if abs(d(cde))<eps1 i=1; j=[1 0];
  while abs(d(cde-i))<eps1 i=i+1; j(i+1)=0; end; d=deconv(d,j);
  cde=cde-i;  end
if abs(n(cnu))<eps1 i=1; j=[1 0];
  while abs(n(cnu-i))<eps1 i=i+1; j(i+1)=0; end; n=deconv(n,j);
  cnu=cnu-i;  end
if d(1)~=1 error('Nenner(1)=1 erwartet!'); end; j=0;
while abs(n(j+1))<eps1 j=j+1; end;   Vz=n(j+1);
for i=(1+j):cnu   n(i)=n(i)/Vz;  end; k=j;
if k==0 
  disp('Warnung !')
  disp('Algorithmus wird auf sprungfaehige Strecke angewendet')
  end
if k>=3 
  disp('Warnung: Ist die Totzeit tatsaechlich so gross ?'); end
for i=(2+j):cnu  siz(1,1+i-j)=n(i); end
siz(1,1)=Vz; siz(1,2)=cnu-1-j;
RZKo=1; DKo=0;
if code(4)==1
  if code(3)==1
    disp('Warnung !')
    disp('Keine optimale Loesung fuer Regelkreis moeglich !')
    end
  if cde==1 
    error('Die eingegebene innere Strecke besitzt keine Pole')
    end
  pz=roots(d); l=cde-1;
  disp('Die Strecke besitzt folgende Pole:'); disp(pz)
  hid=1; shid=1; sRZKo=1;
  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" -sinnvoll da komplexe Pole im Paar kompensiert werden
cde=shid; d=hid; DKo=sRZKo-1; 
end 		%RZKo liegt berechnet vor
if code(4)==2
  if code(3)==1
    disp('Warnung !')
    disp('Keine optimale Loesung fuer Regelkreis moeglich !'); end
  if cde==1 
    error('Die eingegebene innere Strecke besitzt keine Pole'); end
  pz=roots(d); pz=sort(pz); l=cde-1; hid=1; shid=1;
  sRZKo=1; RZKo=1; DKo=0;
  if l>4 
    disp('Warnung: Es werden nur die 4 groessten Pole kompensiert')
    m=l-4; pz1(1:4)=pz(m+1:l); l=4;
    if (imag(pz1(1))~=0)&(imag(pz1(2))==0) pz1(2)=pz(m); end
  else pz1=pz;
  end
  for i=1:l
    sRZKo=sRZKo+1; RZKo(sRZKo)=0; arbvek=1; 
    for j=1:(sRZKo-1)
      arbvek(j+1)=RZKo(j+1)-RZKo(j)*pz1(i);
    end; RZKo=arbvek;
  end; 
cde=shid; d=hid; DKo=sRZKo-1; 
end 		%RZKo liegt berechnet vor
for i=2:cde    siz(2,1+i)=d(i);  end
siz(2,1)=k; siz(2,2)=cde-1;
if ((code(6)==0)&(code(3)==1))
  disp('Warnung !')
  disp('Existenz eines I-Anteiles vor Stoermischstelle ueberpruefen!')
end
if ((code(6)==0)&(siz(2,2)==0))
  error('Ohne Integralanteil (code(6)=0) keine Berechnung moeglich');
  end
if ((code(5)==0)&(code(3)==0)&(DKo>1)&(code(7)==1))
  disp('Warnung !')
  disp('Berechnung mit "gedaempftem" D-Anteil')
  disp('nur in Sonderfaellen brauchbar!!')
end
RZO=DKo+code(5);b1=0;
if RZO>4 disp('Warnung: Ist diese hohe Reglerordnung sinnvoll ?');
end
[R]=havobo(siz,code(3:7));
if code(3)==0       %Fuehrungsoptimierung
  if code(7)==0     %ohne Daempfung des D-Anteiles
    for i=1:(code(5)+1)
      svek(i,1)=-R(i,1);
      for j=1:(code(5)+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 code(5)=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 !')
    disp('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   % Stoeroptimierung
  if code(5)==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
      disp('Keine Berechnung moeglich !')
      error('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])';
    erg=V0(2);       % Wahl der kleineren Verstaerkung(Schwingungen) 
                     % mit "erg=V0(1)" auf gr. Verstaerkung umschaltbar
    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 code(5)==2		% maschinelle Vorgabe der Anfangswerte
    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)=cde-1-j; d(cde)=0;
      for i=3:(cde+1) siz(2,i)=d(i-1); end
    else
      anfv=[siz(1,1) siz(1,1) siz(1,1)]; j=0;
    end
    code(5)=1; [E]=havobo(siz,code(3:7));
    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]); 
      erg=V0(2);     % Wahl der kleineren Verstaerkung(Schwingungen) 
	             % mit "erg=V0(1)" auf gr. Verstaerkung umschaltbar  
      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 %Reglerzaehleranteil code(5) wurde nicht mit 1 oder 2 vorgegeben
    error('Dieser Reglergrad ist nicht vorbereitet')
  end
  rzs=sum(erg);
end
nuco=conv(erg,RZKo);
if code(3)==1  % Filterberechnung
  [j,i1]=size(erg); nufil(i1)=0;  
  nufil(1)=rzs/nuco(1); defil=erg/nuco(1);
else nufil=1; defil=1; end
if ((-code(6)*b1)~=0)     deco=[1 b1-code(6) -code(6)*b1];
elseif ((b1-code(6))~=0)  deco=[1 b1-code(6)];
else                      deco=[1];  end
[j,i1]=size(nuco); [j,i2]=size(deco);
if i1>i2 deco(i1)=0; elseif i1<i2 nuco(i2)=0; end
end
