%
 function [c, ceq, c_Z, ceq_Z]  = con_post(Z,data)
%
% This is the constraint function used by the minimization routine
% 
% Z is the array of current control values
%
% data is a structure that contains information about
% the problem.
%  
% This routine is called after the objective function
% The state and control histories are stored in the
% structure SOLN in the 'base' workspace. 
% This was done after evaluating the trajectory in the
% call to the objective function. Here we need only
% capture the SOLN structure and then evaluate end-point constraints
%
  soln = evalin('base','SOLN');
%
  x_final = soln.x(:,end);
%
  if isfield(data,'Phi_eq')
    ceq   = feval(data.Phi_eq,x_final);
    if nargout > 2
       ceq_Z = grad_nonl(data.Phi_eqj,x_final,data);
    end
  else
    ceq = [];
    ceq_Z = [];
  end
%
  if isfield(data,'Phi_in')
    c   = feval(data.Phi_in,x_final);
    if nargout > 2
       c_Z = grad_nonl(data.Phi_inj,x_final,data);
    end
 else
    c   = [];
    c_Z = [];
  end
%
%
% return
%
%===========================================================
% sub function for con_post
%===========================================================
%
 function C_z = grad_nonl(jacobian,x_final,data)
%
% jacobian is the name of the routine for evaluating the
%          Jacobian of the end-conditions
% x_final is the vector-valued state at the end-time
%
% S is the  triply indexed array
%        S(i,j,k) = \frac{\partial x_N(i)}{\partial u_k(j)} 
%
  S = evalin('base','SENSE'); % retrieve the Sensitivities
%
   nabla_x_C  = feval(jacobian,x_final);
   [N_constr, N_states] = size(nabla_x_C);
%
  tmp = [];
  for k=1:data.N_panels
     tmp =  [ tmp   [nabla_x_C*S(:,:,k)]' ];
   end %  k-loop
%  
  C_z = [reshape(tmp,N_constr,data.N_panels*data.N_controls)]';
%
% return
%
