function [Reg] = dr_d_wedz(rt,Vs,in2,in3,in4)
%function [Reg] = dr_d_wedz(rt,Vs,in2,in3,in4)
%Optimization of PID-type controllers for 5 plant types, based on method:
%Naslin Polynomial Method (Double Ratio) plus pre-filter for delayed inputs
%                  filter_denominator = controller_numerator (ex gain)
%in2,in3,in4: sequence and meaning of inputs --> see plant definition below
%Vs:       plant gain
%rt:       plant and controller type == {1, 2; 3, 4; 5, 6; 7, 8}
%          1:  PID-type general;      plant 1=2: gain, 3 PT1
%           2:    "     double zero;
%          3:  PID-type general;      plant 3=4: gain, 2 PT1, 1 INT
%           4:    "     double zero;
%          5:  PID-type general;      plant 5=6: gain, 1 PT1, 2 INT
%           6:    "     double zero;
%          7:  PID-type general;      plant 7=8: gain, 1 OE,  1 PT1
%           8:    "     double zero;
%          9:  PID-type general;     plant 9=10: gain, 1 OE,  1 INT
%           10:    "    double zero;
%  PT1 - first-order delay-element, parameter: time-constant
%  OE  - oscillating element, parameters:      damping, eigen-time constant
%  INT - integrator element, parameter:        time-constant
%  REG - cell array containing controller parameters,
%        for parameter definitions in cell array REG 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
if ~isreal(Vs)||~isscalar(Vs)  disp('Input Vs improper'); end
if ~isreal(in2)||~isscalar(in2) disp('Input 3 improper'); end
if ~isreal(in3)||~isscalar(in3) disp('Input 4 improper'); end
if ~isreal(in4)||~isscalar(in4) disp('Input 5 improper'); end

%Formulas correspond to manuscript pages 21 till 25
dz=mod(rt,2);
   if ~d_off   %switch off disp commands - useful for nested function calls
if nargout==1 NO=0; else NO=1; end  %for text output only
t1='gain, time-constant1, time-constant2, time-constant3';
t2='gain, time-constant1, time-constant2, integrator time-const1';
t3='gain, time-constant1, integrator time-const1, integrator time-const2,';
t4='gain, damping, eigen-time constant1, time-constant1';
t5='gain, damping, eigen-time constant1, integrator time-constant1';
if   dz==1  disp('your choice: PID-controller without double zeros');
else        disp('your choice: PID-controller with double zeros'); end
jj=round(rt/2); %two controllers each refer to the same plants      
disp(['plant ' num2str(jj) ': ' eval(['t' num2str(jj)]);])
disp(['controller ' 177 ' choice: "-"' 10 09 09 09 ...
      'first line:  Vr  Tn  Tv' 10 09 09 09 'second line: Tr1 Tr2 T0']);
if NO==1 disp(['controller ' 177 ' choice: "+"' 10 09 09 09 ...
      'third line:  Vr  Tn  Tv' 10 09 09 09 'fourth line: Tr1 Tr2 T0']);end
   end         % End hidden mode
%Attention, please note: plant gain Vs is only considered in final step
if     rt==1||rt==2  Tsu=in4;         x=Tsu/in2;            y=Tsu/in3;
                     K2=(1+x^2+y^2);  K3=(1+x+y)^3-8*x*y;   K4=(1+x+y)^4;
       if   dz==1 Tn=8*Tsu*K3/K4;
                  Tr1=4*Tsu*K3/K4*(1+sqrt(-1)*sqrt(2*K2*K4/K3^2-1));
                  Tr2=4*Tsu*K3/K4*(1-sqrt(-1)*sqrt(2*K2*K4/K3^2-1));
       else       Tn=16*Tsu*K2/K3;                           %double zero         
                  Tr1=8*Tsu*K2/K3;                           Tr2=Tr1;
       end
                  Vr=K3/(8*x*y);      Tv=4*Tsu*K2/K3;
                                      T0=2*Tsu*(Tr1/in2)*(Tr2/in3)/K2;

elseif rt==3||rt==4  Tsu=in3;    x=Tsu/in2; y=Tsu/in4;
       if   dz==1 Tn=8*Tsu/(1+x);
                  Tr1=4*Tsu/(1+x)*(1+sqrt(-1)*(1-x)/(1+x));
                  Tr2=4*Tsu/(1+x)*(1-sqrt(-1)*(1-x)/(1+x));
       else       Tn=16*Tsu*(1+x^2)/(1+x)^3;                 %double zero
                  Tr1=8*Tsu*(1+x^2)/(1+x)^3;                 Tr2=Tr1;
       end
                  Vr=(1+x)^3/(8*x*y); Tv=4*Tsu*(1+x^2)/(1+x)^3;
                                      T0=2*Tsu*(Tr1/in2)*(Tr2/in4);

elseif rt==5||rt==6  Tsu=in2;    x=Tsu/in3; y=Tsu/in4;
       if   dz==1 Tn=8*Tsu;
                  Tr1=4*Tsu*(1+sqrt(-1));
                  Tr2=4*Tsu*(1-sqrt(-1));
       else       Tn=16*Tsu;                                 %double zero
                  Tr1=8*Tsu;                                 Tr2=Tr1;
       end
                  Vr=1/(8*x*y);       Tv=4*Tsu;
                                      T0=2*Tsu*(Tr1/in3)*(Tr2/in4);
        
elseif rt==7||rt==8  D=in2;      Tau=in3;   Tsu=in4;    x=Tsu/Tau;
       K2=1+2*(1-2*D^2)*x^2;     K3=(1+2*D*x)^3-8*x^2;  K4=(1+2*D*x)^4;
       if  dz==1  Tn=8*Tsu*K3/K4;
                  Tr1=4*Tsu*(K3/K4)*(1+sqrt(-1)*sqrt(2*K2*K4/K3^2-1));
                  Tr2=4*Tsu*(K3/K4)*(1-sqrt(-1)*sqrt(2*K2*K4/K3^2-1));
       else       Tn=16*Tsu*K2/K3;                          %double zero
                  Tr1=8*Tsu*K2/K3;                           Tr2=Tr1;
       end
                  Vr=K3/(8*x^2);      Tv=4*Tsu*K2/K3;
                                      T0=2*Tsu*(Tr1/Tau)*(Tr2/Tau)/K2;
elseif rt==9||rt==10  D=in2;       Tau=in3;      IZK0=in4;
                      K2=2*D^2-1;  x=IZK0/Tau;
       if  dz==1  Tn=4*Tau/D;
                  Tr1=Tn/2*(1+sqrt(-1)*sqrt(1-1/D^2));
                  Tr2=Tn/2*(1-sqrt(-1)*sqrt(1-1/D^2));
       else       Tn=4*Tau*K2/D^3;                          %double zero
                  Tr1=Tn/2;                                 Tr2=Tr1;
       end
                  Vr=x*D^3;           Tv=Tau*K2/D^3;
                                      T0=(Tr1*Tr2/IZK0)/K2;
else error('dr_d_wedz: Incorrect code for controller type!'),
end                                       %final step:
       Reg={Vr/Vs Tn Tv; Tr1 Tr2 T0*Vs};  %consideration of Vs
end
