function [Reg,Reg2] = bo_d_speso(rt,varargin)
%function [Reg(,Reg2)] = bo_d_speso(rt,varargin)
%Optimization of PID-type controllers for 4 plant types, based on method:
%Magnitude Optimum plus pre-filter (delayed inputs) = Symmetric Optimum, 
%          filter_denominator = controller_numerator (gain factored out)
%varargin: for sequence and meaning of inputs see plant definitions below
%rt:       plant-type = {1, 2, 3, 4}:        plant 1: gain, 2 PT1, 1 INT
%                                            plant 2: gain, 3 PT1
%                                            plant 3: gain, 1 OE,  1 PT1
%                                            plant 4: gain, 1 OE,  1 INT
%  PT1 - first-order delay-element, parameter: time-constant
%  OE  - oscillating element, parameters:      damping, eigen-time constant
%  INT - integrator element, parameter:        time-constant
%  REG - vector containing controller parameters,
%        typically the "useful" solution, because of: choice "-" of "+/-"
%  REG2- typically the "unfit"  solution, because of: choice "+" of "+/-"
%        for parameter definitions in vectors REG/REG2 see display command,
%        >>>thus do not deactivate display command for first use<<<
%                      parameter meaning: GC(p)=Vr*[1+1/(p*Tn)+p*Tv]
%                      VR=TN/TR0, p: Laplace variable, thus alternatively:
%                      GC(p)=[1+p*Tn+p^2*Tv*Tn]/(p*T0) OR AGAIN
%                      GC(p)=[1+p*Tr1]*[1+p*Tr2]/[p*T0]
%        hence pre-filter
%               GFi(p)=1/[1+p*Tn+p^2*Tv*Tn]=1/[1+p*(Tr1+Tr2)+p^2*(Tr1*Tr2)]
%  For deactivation of display commands use hidden mode "rt+10" as type
%                                        - useful for nested function calls

%Check of inputs:
if round(rt)~=rt|~isscalar(rt) disp('Input rt improper'); end
if rt>10 rt=rt-10; d_off=1; else d_off=0; end  %realizes "disp" hidden mode
%Formulas correspond to book chapter10.4
varein=varargin;
for j=1:4 if ~isreal(varein{j})||~isscalar(varein{j})
                        disp(['Input ' num2str(j+1) ' improper']); end, end
if rt>0 && rt<5 if nargout==1 NO=0; else NO=1; end %to output unfit result?
                [Reg,Reg2]=PID_d(rt,NO,d_off,varein);
else            error('bo_d_speso: Incorrect code for controller type!'),
end

function [Reg,Reg2]=PID_d(rt,NO,d_off,varein)
   if ~d_off   %switch off disp commands - useful for nested function calls
t1='time-constant 1, time-constant 2, integrator time-constant';       
t2='time-constant 1, time-constant 2, time-constant 3';
t3='damping, eigen-time constant, time-constant';
t4='damping, eigen-time constant, integrator time-constant';
disp('your choice: PID-controller for delayed inputs');
disp(['plant ' num2str(rt) ': gain, ' eval(['t' num2str(rt)]);])
disp(['controller ' ...
      ' first line (' 177 ' choice: "-"):  Vr  Tn  T0  Tv' 10 09 09 09 ...
      'second line (' 177 ' choice: "+"): Vr  Tn  T0  Tv']);
if NO==1 disp([ 09 09 09 ...
      'third line (' 177 ' choice: "-"):  Tr1 Tr2' 10 09 09 09 ...
      'fourth line (' 177 ' choice: "+"): Tr1 Tr2']);end
   end
% case dependen definitions
if rt==1
Vs =varein{1};  T0=varein{4};   Ts2=varein{2};     Tsu=varein{3};
x=Tsu/T0;       y=Tsu/Ts2;      xq=x*y;  
Krt=1+y;        Kq=1+y^2;       Tf0=Tsu;           Ks=0;
elseif rt==2
Vs =varein{1};  Ts1=varein{2};  Ts2=varein{3};     Tsu=varein{4};
x=Tsu/Ts1;      y=Tsu/Ts2;      xq=x*y;  
Krt=1+x+y;      Kq=1+x^2+y^2;   Tf0=Tsu;           Ks=4;
elseif rt==3 
Vs =varein{1};  Dam=varein{2};  Tau=varein{3};     Tsu=varein{4};
x=Tsu/Tau;                      xq=x^2;  
Krt=1+2*Dam*x;  Kq=1+2*xq*(2*Dam^2-1);             Tf0=Tsu;     Ks=4;
elseif rt==4
Vs =varein{1};  Dam=varein{2};  Tau=varein{3};     T0=varein{4};
x=T0/Tau;                       xq=x^2;  
Krt=2*Dam*x;    Kq=2*xq*(2*Dam^2-1);               Tf0=T0;      Ks=0;
end

% equations valid for all cases, except for yq_1/2
FP=(2+sqrt(2));                    FM=(2-sqrt(2));
Vr1=(FP*Krt^3-Ks*xq)/(4*xq);       Vr2=(FM*Krt^3-Ks*xq)/(4*xq);
Tn1=4*Tf0*(FP*Krt^3-Ks*xq)/(FP^2*Krt^4);
                                   Tn2=4*Tf0*(FM*Krt^3-Ks*xq)/(FM^2*Krt^4);
Tv1=2*Tf0*Kq/(FP*Krt^3-Ks*xq);     Tv2=2*Tf0*Kq/(FM*Krt^3-Ks*xq);
Re1=2*Tf0*(FP*Krt^3-Ks*xq)/(FP^2*Krt^4);
                                   Re2=2*Tf0*(FM*Krt^3-Ks*xq)/(FM^2*Krt^4);
Im1=sqrt((2*FP^2*Krt^4*Kq)/(FP*Krt^3-Ks*xq)^2-1);
                          Im2=sqrt((2*FM^2*Krt^4*Kq)/(FM*Krt^3-Ks*xq)^2-1);
Tr1a=Re1*(1+i*Im1);                Tr2a=Re1*(1-i*Im1);
Tr1b=Re2*(1+i*Im2);                Tr2b=Re2*(1-i*Im2);
   if rt<3
      if   rt==2 yq_1=(Tr1a*Tr2a)/(Ts1*Ts2); yq_2=(Tr1b*Tr2b)/(Ts1*Ts2);
      else       yq_1=(Tr1a*Tr2a)/(T0*Ts2);  yq_2=(Tr1b*Tr2b)/(T0*Ts2); end
   else
      yq_1=(Tr1a*Tr2a)/Tau^2;      yq_2=(Tr1b*Tr2b)/Tau^2;     end
T0a=2*Tf0*yq_1/Kq;                 T0b=2*Tf0*yq_2/Kq;
% arrange results
Reg =[ Vr2/Vs  Tn2  T0b*Vs  Tv2; Vr1/Vs  Tn1  T0a*Vs  Tv1 ];
Reg2=[ Tr1b    Tr2b;             Tr1a    Tr2a             ];
