function [Gamma,gamma] = ppp_output_constraint (A,B,C,D,x_0,A_u,Tau,Min,Max,Order,i_y)

  ## usage:  [Gamma,gamma] = ppp_output_constraint (A,B,C,D,x_0,A_u,Tau,Min,Max,Order)
  ##
  ## Derives the output constraint matrices Gamma and gamma
  ## For Constraints Min and Max at times Tau
  ## Initial state x_0 
  ## Order=0 - output constraints
  ## Order=1 - output derivative constraints
  ## etc
  ## NOTE You can stack up Gamma and gamma matrices for create multi-output constraints.

  ## Copyright (C) 1999 by Peter J. Gawthrop
  ## 	$Id$	

  ## Sizes
  [n_x,n_u,n_y] = abcddim(A,B,C,D); # System dimensions
  [n_U,m_U] = size(A_u);	# Number of basis functions
  [n,n_tau] = size(Tau);		# Number of constraint times
  
  if n_tau==0			# Nothing to be done
    Gamma = [];
    gamma = [];
    return
  endif

  ## Defaults
  if nargin<10
    Order = zeros(1,n_tau);
  endif

  if nargin<11
    i_y = 1;			# First output
  endif

  if n != 1
    error("Tau must be a row vector");
  endif
  
  n = length(Min);
  m = length(Max);
  o = length(Order);

  if (n != n_tau)||(m != n_tau)||(o != n_tau)
    error("Tau, Min, Max and Order must be the same length");
  endif
  

  ## Compute Gamma 
  Gamma = [];
  zero_x = zeros(size(x_0));
  for i=1:n_U
    U = zeros(n_U,1); U(i,1) = 1; # Set up U_i
    y_i = ppp_ystar (A,B,C,D,zero_x,A_u,U,Tau); # Compute y* for ith input for each tau
    y_i = y_i(i_y,:); # Pluck out output i_y
    ## Gamma = [Gamma [-y_i';y_i']]; # Put in parts for Min and max
    Gamma_i = [];
    if (Min>-inf)
      Gamma_i = [Gamma_i; -y_i']; # Min part of column
    endif
    if (Max<inf)
      Gamma_i = [Gamma_i;  y_i'];	# Add column for Max
    endif
    Gamma = [Gamma Gamma_i]; # Put in parts for Min and max
  endfor

  ## Compute gamma
  zero_U = zeros(size(U));
  y_x = ppp_ystar (A,B,C,D,x_0,A_u,zero_U,Tau);	# Output if U is zero
  gamma = [];
  if (Min>-inf)
    gamma = [gamma; -(Min-y_x)'];
  endif
  if (Max<inf)
    gamma = [gamma; (Max-y_x)'];
  endif
  
endfunction