function  ppp_1(Name,Inputs,Modes);

# Name = "CantileverBeam"
# Inputs = 1:2
# Modes = [1 2;1 3]

  path(path,"~/Research/CGPC/PPP");  
  
  Name,Inputs,Modes

  ## System
  system(sprintf("mtt -q %s numpar m", Name));
  system(sprintf("mtt -q %s sm m", Name));
  par = eval(sprintf("%s_numpar;",Name));
  eval(sprintf("[A,B,C,D]=%s_sm(par);",Name));
  [n_x,n_u,n_y] = abcddim(A,B,C,D)


  ## Change B
  B = B(:,Inputs);
  [junk,n_u]=size(B);
  n_u

  ## Redo C and D to reveal ALL velocities
   n_y = n_x/2
   C = zeros(n_y,n_x);
   for i = 1:n_y
     C(i,2*i-1) = 1;
   endfor

  ## Sort out D
  D = zeros(n_y,n_u);

  e = eig(A);			# Eigenvalues
  N = length(e);
  frequencies = sort(imag(e));
  frequencies = frequencies(N/2+1:N); # Modal frequencies

  ## Controller
  ## Controller design parameters
  t = [0.9:0.01:1.0];		# Optimisation horizon

  ## Specify input basis functions 
  ##  - damped sinusoids with same frequencies as beam
  damping_ratio = [0.1 0.1];		# Damping ratio of inputs
  A_u=[];
#   Modes = [1  3  
# 	   1  2];		# Choose modes to be controlled by each input 
  if n_u == 1			# Put all modes on each input
    Modes = [Modes; Modes];
  endif
  
  Modes

  for i=Inputs
    A_ui = damped_matrix(frequencies(Modes(i,:)),damping_ratio(i)*ones(size(Modes(i,:))));
    A_u = [A_u;A_ui];
    u_0 = ones(n_x,1);		# Initial conditions
  endfor
  

  A_w = 0;			# Setpoint
  Q = ones(n_y,1);		# Output weighting

  ## Design 
  disp("Control design - unconstrained");
  [k_x,k_w,K_x,K_w,Us0,J_uu,J_ux,J_uw] = ppp_lin  (A,B,C,D,A_u,A_w,t); # Unconstrained design

  ## Organise some plots
   dt = 0.002;		# Time increment
#   T =  0:dt:t(length(t));	# Time starting at zero but past horizon
   T =  0:dt:1.0;		# Time 

  ## Set up an "typical" initial condition
  x_0 = zeros(n_x,1);
  x_0(2:2:n_x) = ones(1,n_x/2);	# Set initial twist to 1.
#  disp("Computing OL response");
#  [Ys,Us] = ppp_transient (T,A_u,-(K_x*x_0)',u_0,A,B,C,D,x_0); # Compute open-loop control
  
  
  disp("Computing OL response (no control)");
  [Y_0] = ppp_sm2sr(A, B, C, D, T, zeros(n_u,1), x_0); # Compute Closed-loop control
  title("Y (no control)");
  xlabel("Time")
  grid; 
  plot(T,Y_0);
  psfig(sprintf("%s_pppy0",Name));

  disp("Computing Unconstrained CL response");
  [y_u,X_u] = ppp_sm2sr(A-B*k_x, B, C, D, T, zeros(n_u,1), x_0);
  u_u = -k_x*X_u';


  title("Unconstrained closed-loop control - y");
  xlabel("Time")
  grid; 
  plot(T,y_u);
  psfig(sprintf("%s_pppy",Name));

  title("Unconstrained closed-loop control - u");
  xlabel("Time")
  grid; 
  plot(T,u_u);
  psfig(sprintf("%s_pppu",Name));

#   ## Constrained version
  
#   ## Constraints - u
#   Tau_u = [0:0.01:1];
#   one = ones(size(Tau_u));
#   limit = 1e10;
#   Min_u = -limit*one;
#   Max_u =  limit*one;
#   Order_u = 0*one;

#   ## Constraints - y
#   Tau_y = [];
#   one = ones(size(Tau_y));
#   limit = 1e5; 
#   Min_y = -limit*one;
#   Max_y =  limit*one;
#   Order_y = 0*one;

#   ## Simulation

#   ## Constrained - open-loop
#   Gamma = [];
#   gamma = [];
#   for i=1:n_u			# Put same constraints on each input
#     [Gamma_i,gamma_i] = ppp_input_constraint (A_u,Tau_u,Min_u,Max_u,Order_u,i,n_u);
#     Gamma = [Gamma; Gamma_i];
#     gamma = [gamma; gamma_i];
#   endfor


#   disp("Open-loop simulations...");
#   ## Constrained OL simulation
#   W = zeros(n_y,1);
#   [u,U] = ppp_qp (x_0,W,J_uu,J_ux,J_uw,Us0,Gamma,gamma);
#   T = [0:t(2)-t(1):t(length(t))];
#   [ys,us] = ppp_ystar (A,B,C,D,x_0,A_u,U,T);

#   title("Constrained y*");
#   xlabel("t");
#   grid;
#   plot(T,ys)

#   ## Unconstrained OL simulation
#   [uu,Uu] = ppp_qp (x_0,W,J_uu,J_ux,J_uw,Us0,[],[]); 
#   [ysu,usu] = ppp_ystar (A,B,C,D,x_0,A_u,Uu,T);

#   title("Unconstrained y*");
#   xlabel("t");
#   grid;
#   plot(T,ysu)

#   ## Non-linear - closed-loop
#   disp("Closed-loop simulation");
#   Tc = [0:4e-4:1.0];	
#   Delta_ol = 1.0;
#   [yc,uc,J] = ppp_qp_sim (A,B,C,D,A_u,A_w,t,Q, T, Tau_u,Min_u,Max_u,Order_u, Tau_y,Min_y,Max_y,Order_y,W,x_0,Delta_ol);

#   title("y,y*,u and u*");
#   xlabel("t");
#   grid;
#   plot(T1,y,T,ys,T1,u,T,us);

endfunction