function [ssc]=StateSpaceControlDesign(mot,mode)
  % !!! first it need to run: [mot1,mot2]=identifyFxFyStage() to build a motor object !!!
  %
  % builds a state space controller designed for the plant.
  % shows step answers of open and closed loop, also for the observer controller and the final discrete observer
  %
  % the matchich simulink model is: 'observer'

  %References:
  %http://ctms.engin.umich.edu/CTMS/index.php?example=Introduction&section=ControlStateSpace
  %space state controller:
  % web(fullfile(docroot, 'simulink/examples.html'))
  % web(fullfile(docroot, 'simulink/examples/inverted-pendulum-with-animation.html'))
  % web(fullfile(docroot, 'simulink/examples/double-spring-mass-system.html'))
  %
  % https://www.youtube.com/watch?v=Lax3etc837U

  %plant and model
  % ss_plt :real plant (model of real plant)
  % ss_c1  :current, mechanic, no resonance
  % ss_d1  :simpl. current, mechanic, no resonance
  % ss_1  :no current, mechanic, no resonance
  % ss_0  :no current, simpl. mechanic, no resonance

  %prefilt:prefilter mode
  %0 no filter
  %1 inverse resonance filter
  %2 manual setup filter

  %verb: mode(bits) to plot/simulate
  % 0:  1: poles of model and placed poles of controller
  % 1:  2: bode plots of open loop
  % 2:  4: step answer on open loop
  % 3:  8: step answer on closed loop with space state controller
  % 4: 16: step answer on closed loop with observer controller
  % 5: 32: step answer on closed loop with disctrete observer controller
  % 6: 64: plot all closed loop bode and pole-zero diagrams of desPos->actPos
  % 7:128: bode plot of filt_pos_err


  %use_lqr: use lqr instead of pole placement
  t=0:1E-4:.5; %time vector for simulations
  verb=0;
  use_lqr=0;
  MaxDac=2011.968;
  tfDesPos=tf(1,1);
  tfPosErr=tf(1,1);
  ssc=struct();
  ssc.sel1={3,[3]};
  %locate poles: 2500rad/s = 397Hz, 6300rad/s = 1027Hz
  switch mode
    case -1 %TESTING
      
    case 0 %for document: best observer without prefilter
      ss_plt=mot.ss_plt; %ss_plt ss_c1 ss_d1 ss_1 ss_0
      ss_mdl=mot.ss_cp;  %ss_plt ss_c1 ss_d1 ss_1 ss_0
      if mot.id==1
        pl=[-3300 -3200 -2900 -2800];
      else
        pl=[-3300 -3200 -2700 -2600];
      end
      plObs=2*pl;
    case 1 %for document: best observer without prefilter
      ss_plt=mot.ss_plt; %ss_plt ss_c1 ss_d1 ss_1 ss_0
      ss_mdl=mot.ss_cp;  %ss_plt ss_c1 ss_d1 ss_1 ss_0
      if mot.id==1
        pl=[-3300 -3200 -2900 -2800];
      else
        pl=[-3300 -3200 -2700 -2600];
      end
      plObs=2*pl;
      tfDesPos=Prefilt(mot,2);%user designed envelope
      tfPosErr=Prefilt(mot,1);%inverse resonance
    case 2 %for document: best observer with prefilter
      ss_plt=mot.ss_plt; %ss_plt ss_c1 ss_d1 ss_1 ss_0
      ss_mdl=mot.ss_dp;  %ss_plt ss_c1 ss_d1 ss_1 ss_0
      if mot.id==1
        pl=[-2200 -2100 -2000]; % stable with scaling of .05 .. 1.0;
      else
        pl=[-2500 -900 -800];
      end
      plObs=2*pl;
    case 3
      ss_plt=mot.ss_plt; %ss_plt ss_c1 ss_d1 ss_1 ss_0
      ss_mdl=mot.ss_dp;  %ss_plt ss_c1 ss_d1 ss_1 ss_0
      if mot.id==1
        pl=[-2200 -2100 -2000]; % stable with scaling of .05 .. 1.0;
      else
        pl=[-2500 -900 -800];
      end
      plObs=2*pl;
    case 4
      %this is the hovering ball model
      A = [ 0   1   0; 980  0  -2.8;0   0  -100 ];
      B = [ 0 0 100 ]';
      C = [ 1 0 0 ];
      D = 0;
      ssBall=ss(A,B,C,D,'InputName','iCmd','OutputName','actPos');
      ss_plt=ssBall; %ss_plt ss_c1 ss_d1 ss_1 ss_0
      ss_mdl=ssBall;  %ss_plt ss_c1 ss_d1 ss_1 ss_0
      pl=[-10+10i -10-10i -50];
      plObs=[-100 -101 -102];
  end

  [Am,Bm,Cm,Dm]=ssdata(ss_mdl);
  
  if bitand(verb,1) && use_lqr==0
    format compact
    format shortG
    disp(pole(ss_mdl)) %==eig(Am)
    %damp(ss_mdl) %further informations
    disp(pl)
    disp(plObs)
    format short
  end  
  
  if bitand(verb,2)
    figure();h=bodeplot(ss_plt,ss_mdl);
    setoptions(h,'IOGrouping','all')
  end

  xp0 = zeros(1,length(ss_plt.A));
  xm0 = zeros(1,length(Am));

  if bitand(verb,4)
    % step answer on open loop:
    u = ones(size(t));
    [yp,t,x] = lsim(ss_plt,u,t,xp0);
    [ym,t,x] = lsim(ss_mdl,u,t,xm0);
    figure();hold on; plot(t,yp,'DisplayName',ss_plt.OutputName{1}) 
    plot(t,ym,'--','DisplayName',ss_plt.OutputName{1});
    title('step on open loop (plant and model)');
    %legend('plt.iqMeas','plt.iqVolts','plt.actPos','mdl.iqMeas','mdl.iqVolts','mdl.actPos')
    legend('location','best')
  end
  %w0=abs(poles);
  %ang=angle(-poles);
  %-------------------
  %p=w0.*exp(j.*ang)

  % *** space state controller ***
  %
  %place poles for the controller feedback
  if use_lqr %use the lqr controller
    Q=eye(length(Am));
    R=1;
    [K,P,E]=lqr(Am,Bm,Q,R,0);
  else
    K = place(Am,Bm,pl);
    %K = acker(Am,Bm,pl);
  end %if lqr

  V=-1./(Cm*(Am-Bm*K)^-1*Bm); %(from Lineare Regelsysteme2 (Glattfelder) page:173 )
  if length(V)>1
    disp(['scaling (should be actPos): ' ss_mdl.OutputName{end}])
    V=V(end); % only the position scaling needed
  end

  ss_cl = ss(Am-Bm*K,Bm*V,Cm,0,'Name','space state controller','InputName',ss_mdl.InputName,'OutputName',ss_mdl.OutputName);
  if bitand(verb,8)
    % step answer on closed loop with space state controller:
    [y,t,x]=lsim(ss_cl,V*u,t,xm0);
    figure();plot(t,y);title('step on closed loop');
  end

  % *** observer controller ***
  %
  %observer poles
  if use_lqr %use the lqr controller
    Q=eye(length(Am'));   % ??????????????? CHANGES NEEDED ????????????
    R=eye(size(Cm,1));
    [L,P,E]=lqr(Am',Cm',Q,R,0);
  else
    L=place(Am',Cm',plObs)';
    %L=acker(A',C',plObs)';
  end

  At = [ Am-Bm*K             Bm*K
         zeros(size(Am))    Am-L*Cm ];
  Bt = [    Bm*V
         zeros(size(Bm)) ];
  Ct = [ Cm    zeros(size(Cm)) ];

  Dt=0;
  ss_t = ss(At,Bt,Ct,Dt,'Name','observer controller','InputName',{'desPos'},'OutputName',ss_mdl.OutputName);
  if bitand(verb,16)
    % step answer on closed loop with observer controller:
    figure();lsim(ss_t,ones(size(t)),t,[xm0 xm0]);title('step on closed loop with observer');
  end

  % *** disctrete observer controller ***
  %
  Ts=1/5000; % 5kHz
  ss_tz = c2d(ss_t,Ts);
  [Atz,Btz,Ctz,Dtz]=ssdata(ss_tz );
  ss_tz.Name='discrete obsvr ctrl';

  if bitand(verb,32)
    % step answer on closed loop with disctrete observer controller:
    t = 0:Ts:.05;
    figure();lsim(ss_tz ,ones(size(t)),t,[xm0 xm0]);title('step on closed loop with observer discrete');
  end

  if bitand(verb,64)
    %plot all bode diagrams of desPos->actPos
    figure();
    idx=length(ss_cl.OutputName);
    h=bodeplot(ss_cl(idx),ss_t(idx),ss_tz(idx));
    setoptions(h,'FreqUnits','Hz','Grid','on');legend('location','sw');

    figure();
    h=pzplot(ss_cl(idx),ss_t(idx));
    setoptions(h,'FreqUnits','Hz','Grid','on');legend('location','sw');
    figure();
    h=pzplot(ss_tz(idx));
    setoptions(h,'FreqUnits','Hz','Grid','on');legend('location','sw');
  end

  %calculate matrices for the simulink system
  Ao=Am-L*Cm;
  Bo=[Bm L];
  Co=K;
  Do=zeros(size(Co,1),size(Bo,2));

  ss_o = ss(Ao,Bo,Co,Do,'Name','observer controller','InputName',[{'desPos'}; ss_mdl.OutputName ],'OutputName',{'k*xt'});

  %discrete plant
  %ss_pltz = c2d(ss_plt,Ts);
  %[Apz,Bpz,Cpz,Dpz]=ssdata(ss_pltz);

  %discrete observer controller
  ss_oz = c2d(ss_o,Ts);

  %discrete prefilter
  tfDesPos_z=c2d(tfDesPos,Ts);
  tfPosErr_z=c2d(tfPosErr,Ts);

  if bitand(verb,128)
    h=bodeplot(filt_pos_err,filt_pos_err_z);
    setoptions(h,'FreqUnits','Hz','Grid','on');legend('location','sw');
  end

  %state space controller
  for k=["Ts","ss_plt","ss_o","ss_oz","tfDesPos","tfDesPos_z","tfPosErr","tfPosErr_z","V","MaxDac"]
    ssc=setfield(ssc,k,eval(k));
  end

  mat2py=struct();
  %[ozA,ozB,ozC,ozD]=ssdata(ss_oz);
  %[pos_err_num,pos_err_den]=tfdata(filt_pos_err_z);

  mat2py.Ts=Ts;
  mat2py.V=V;
  mat2py.MaxDac=MaxDac;
  mat2py.ozA=ss_oz.A;
  mat2py.ozB=ss_oz.B;
  mat2py.ozC=ss_oz.C;
  mat2py.ozD=ss_oz.D;
  mat2py.ozInpName=ss_oz.InputName;
  mat2py.ozOutName=ss_oz.OutputName;
 
  fn=[pwd '/' sprintf( 'ssc%d.mat',mot.id)];
  save(fn,'-struct','mat2py');
  disp(['saved ' fn]);
end

function pf=Prefilt(mot,mode)
  switch mode
    case 0 %no filter
        pf=tf(1,1);
    case 1 %inverse resonance
      if mot.id==1
        den=mot.mdl.num1;%num=1;
        num=mot.mdl.den1;%den=[1 0 0];
        pf=tf(num,den);
      else
        den=conv(conv(conv(mot.mdl.num1,mot.mdl.num2),mot.mdl.num3),mot.mdl.num4);%num=1;
        num=conv(conv(conv(mot.mdl.den1,mot.mdl.den2),mot.mdl.den3),mot.mdl.den4);%den=[1 0 0];
        pf=tf(num,den);
      end
    case 2
      if mot.id==1
        %f=200;w0=f*2*pi;  num1=[1 300 w0^2]; den1=[1 200 w0^2];
        %numV=num1;
        %denV=den1;
        %pf=tf(numV,denV);
        %Lag
        f=[200 300]; w=f*2*pi; T=1./w;
        pf=tf([T(1) 1],[T(2) 1]);
      else
        %Lag
        f=[200 400]; w=f*2*pi; T=1./w;
        tf1=tf([T(1) 1],[T(2) 1]);
        %bo = bodeoptions;
        %bo.FreqUnits = 'Hz'; bo.MagUnits='abs'; bo.Grid='on';
        %bode(tf1,bo)

        %k=1.2; aa=2; f=[40 60];w=f*2*pi;   tf([1 33 w0^2];  den3=[1 20 w0^2];
        %f=277;w0=f*2*pi;  num1=[1 20 w0^2];  den1=[1 500 w0^2];
        %f=138;w0=f*2*pi;  num2=[1 300 w0^2]; den2=[1 100 w0^2];
        %f=60;w0=f*2*pi;   num3=[1 33 w0^2];  den3=[1 20 w0^2];
        %numV=conv(conv(num1,num2),num3);
        %denV=conv(conv(den1,den2),den3) ;
        %pf=tf(numV,denV);
        pf=tf1;
      end
  end
  %controlSystemDesigner('bode',1,pf);    % <<<<<<<<< This opens a transferfunction that can be edited
end