%
 function [ val, val_Z] = obj_post(Z,data)
%
% This is the objective function used by the minimization routine
% 
% Z is the array of current control values
%
% data is a structure that contains information about
% the problem.
%  
  soln.u            = reshape(Z, data.N_controls, data.N_panels);
  [ soln.t, soln.x] = generate(data,soln.u);
%
  val    = feval(data.Mayer,soln.x(:,1+data.N_panels));
%
% For later use we set the SOLN structure in the base workspace
%
  assignin('base','SOLN',soln);
%
% check to see of gradient is also required
%
  if ( nargout > 1 )
    val_Z = grad_cost(soln.x(:,end),data);
  end
%
% return
%
%========================================================
% internal subfunctions for obj_post
%========================================================
 function [t_array, x_array] = generate(data,u_array)
%
% This is the function used to generate the
% values of the vector states at the N_panels + 1 points
% starting at t = 0, to t = T_f
% 
% u_array is the array of current control values
%
% data is a structure that contains information about
% the problem.
%  
  x_array = zeros(data.N_states,1+data.N_panels);
  x_array(:,1) = data.X_initial(:);
  t_array      = (data.T_f / data.N_panels)*[0:1:data.N_panels];
%
 for k=1:data.N_panels
    [T, Y]         = ode23(data.Rhs,[t_array(k) t_array(k+1)], ...       
                      x_array(:,k)',[],u_array(:,k));
    x_array(:,k+1) = Y(end,:)'; % assign the last row
 end
%
% return
%
%========================================================
%
 function F_z = grad_cost(x_final,data)
%
% data.Mayerj is the name of the routine for evaluating the
%        gradient of the Mayer functional
%
% 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)} 
%
  soln = evalin('base','SOLN');
  S    = sense(soln,data);
%
  nabla_x_F = feval(data.Mayerj,x_final);
  nabla_x_F = nabla_x_F(:)'; % Make sure it's a row
%
  tmp = zeros(data.N_controls,data.N_panels);
  for k=1:data.N_panels
     tmp(:,k) =  [nabla_x_F*S(:,:,k)]';
   end %  k-loop
%  
  F_z = reshape(tmp,1,data.N_panels*data.N_controls);
%
 assignin('base','SENSE',S); % save the Sensitivities
%
% return
%
%========================================================
%
 function S = sense(soln,data)
%
% The purpose of this routine is to generate the sensitivity maps
% (or partial derivatives) from the control (values) to the state values. 
% The optimization problem requires only the sensitivity of the final values of the states.
% We have discretized the control function to be piecewise constant. 
% In the following we refer to the  interval  [t_{i-1} t_i] as the i-th grid panel. 
% On this panel we have the control vector constant at the value u_i. 
% The required sensitivity maps are 'built-up' from local maps on each panel.
%
% The semi-group nature of the solution of the differential equation implies that
% the control on the j-th panel can have no effect on state-values at previous panels.
%  Suppose at a fixed j-value we have the sensitivities for 
% \frac{\partial x_{j-1}}{\partial u_k} \equiv  [ S^j_k ]   k = 1, ..., j-1.
% Note, these are  N_states by N_controls matrices.
% We determine the 'single-panel' maps
% \frac{\partial x_{j}}{\partial u_j} \equiv  [ \Gamma^j  ]
% (an N_states by N_controls matrix)
% and 
% \frac{\partial x_{j}}{\partial x_{j-1}} \equiv [ \Phi^j ]
% (an N_states by N_states matrix).
% We can then generate the collection of sensitivities
% [ S^{j+1}_k ]  for k = 1, ..., j, from the rules
% [ S^{j+1}_k ]  = \Phi^j [ S^j_k ], k = 1, ..., j-1, and
% [ S^{j+1}_j ]  = \Gamma^j .
% We move across the panels from left to right and end with
%  [ S^{N}_j ]  =  \frac{\partial x_{N}}{\partial u_j}, j = 1,..N
%  (these are  N_states by N_controls matrices).
%
% The S-matrices are stored in the triply indexed array
%  S(i,j,k) 
% i = 1:N_states
% j = 1:N_controls
% k = 1:N_panels
%
% pre-allocate storage for the array
  S = zeros(data.N_states,data.N_controls,data.N_panels);
%
  for k=1:data.N_panels
%
% integrate for local sensitivities
% for this panel
%
      local.t_l   = soln.t(k);
      local.x_l   = soln.x(:,k);
      local.t_r   = soln.t(k+1);
      local.x_r   = soln.x(:,k+1);
      local.u     = soln.u(:,k);
      PG0         = [eye(data.N_states) ...
                     zeros(data.N_states,data.N_controls)];
      [ T, PG ] = ode23('rhs_post',[soln.t(k), soln.t(k+1)], ...
                         PG0,[],local,data);
      tmp         = reshape(PG(end,:),data.N_states, ...
                             data.N_states+data.N_controls);
%
%     Phi         = tmp(:,1:data.N_states);
%     Gam         = tmp(:,1+data.N_states:end);
%
%  propagate sensitivities across the panel
%
      S(:,:,k) = tmp(:,1+data.N_states:end); % Gam 
	  l = 1;
	  while l < k
	     S(:,:,l)=tmp(:,1:data.N_states)*S(:,:,l);% = Phi*S(:,:,l)
	     l = l + 1;	
	  end %   l-loop
%
  end % k-loop
