exec('/store3/sci/scilab-2.2/demos/control/scheme.dem');
s=poly(0,'s');z=poly(0,'z');
x_message('Example of PID Design ')
n=x_choose(['Continuous time';'Discrete time'],'Select time domain');
select n
 case 0
  warning('Demo stops!');return;
 case 1
  dom='c';
  s=poly(0,'s');
  str='[(s-1)/(s^2+5*s+1)]';
  rep=x_dialog('Nominal plant?',str)
  if rep==[] then return,end
  Plant=evstr(rep); 
  Plant=syslin('c',Plant);
 case 2
  dom='d'
  z=poly(0,'z');
  str='(z+1)/(z^2-5*z+2)'
  rep=x_dialog('Nominal plant?',str)
  if rep==[] then return,end
  Plant=evstr(rep)
  Plant=syslin('d',Plant);
end   
 //Nominal Plant
P22=tf2ss(Plant);    //...in state-space form
[ny,nu,nx]=size(P22);
   defv=['-1.2','1','0.1'];
   if dom='d' then defv=['-10','1','0.1'];end
while %t
   if dom='c' then
   title='Enter your PID controller K(s)=Kp*(1+T0/s+T1*s)';
   end
   if dom='d' then
   title='Enter your PID controller K(z)=Kp*(1+T0/z+T1*z)';
   end
   defv=x_mdialog(title,['Kp=';'T0=';'T1='],defv);
   if defv=[] then warning('Demo stops!');return;end
   Kp=evstr(defv(1));T0=evstr(defv(2));T1=evstr(defv(3));
   if dom='c' then
   Kpid=tf2ss(Kp*(1+T0/s+T1*s));
   end
   if dom='d' then
   Kpid=tf2ss(Kp*(1+T0/z+T1*z));
   end
   W=[1, -P22;
      Kpid,1];Winv=inv(W);

   disp(spec(Winv(2)),'closed loop eigenvalues');//Check internal stability
   if maxi(real(spec(Winv(2)))) > 0 then 
     x_message('You loose: closed-loop is UNSTABLE!!!');
                                    else
     x_message('Congratulations: closed-loop is STABLE !!!');
     break;
   end
end

[Spid,Rpid,Tpid]=sensi(P22,Kpid);  //Sensitivity functions
Tpid(5)=clean(Tpid(5));

disp(clean(ss2tf(Spid)),'Sensitivity function');
disp(clean(ss2tf(Tpid)),'Complementary sensitivity function');

resp=['Frequency response';'Time response'];
while %t do
n=x_choose(resp,'Select response(s)');
if degree(Tpid(5))>0 then
    error('Improper transfer function!')
end
Tpid(5)=coeff(Tpid(5));
select n
 case 0
  break
 case 1
  xbasc(1);xset("window",1);xselect();bode(Tpid);
 case 2
  if Plant(4)='c' then
   defv=['0.1','50'];
   title='Enter Sampling period and Tmax';
   rep=x_mdialog(title,['Sampling period?';'Tmax?'],defv)
   if rep==[] then break,end
   dttmax=evstr(rep);
   dt=evstr(dttmax(1));tmax=evstr(dttmax(2));
   t=0:dt/5:tmax;
   n1=x_choose(['Step response?';'Impulse response?'],'Simulation:');
   if n1==0 then
      warning('Demo stops!');return;
   end
   if n1==1 then 
      xbasc(1);xset("window",1);xselect();
      plot2d([t',t'],[(csim('step',t,Tpid))',ones(t')])
   end
   if n1==2 then
     xbasc(1);xset("window",1);xselect();
     plot2d([t',t'],[(csim('impul',t,Tpid))',0*t'])
   end
  end
  if Plant(4)='d' then
   defv=['100'];
   title='Tmax?'
   rep=x_mdialog(title,['Tmax='],defv)
   if rep==[] then break,end
   tmax=evstr(rep);
   while %t do
   n=x_choose(['Step response?';'Impulse response?'],'Simulation:');
   select n
   case 0 then
     break
   case 1 then
      u=ones(1,Tmax);u(1)=0;
      xbasc(1);xset("window",1);xselect();
      plot2d([(1:tmax)',(1:tmax)'],[(dsimul(Tpid,u))',(ones(1:tmax)')])
   case 2 then
      u=zeros(1,Tmax);u(1)=1;
      xbasc(1);xset("window",1);xselect();
      plot2d((1:Tmax)',(dsimul(Tpid,u))')
   end
   end
  end
end
end