function test01 ( )

%*****************************************************************************80
%
%% test01() iterates once with two trial controls for a model problem.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    01 September 2022
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test01():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Test optimization_biological()\n' );

  for ufun = 1 : 2

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Using control function U #%d\n', ufun );
  fprintf ( 1, '\n' );
%
%  Integrate forward in time for X.
%
    t0 = 0.0;
    t1 = 5.0;
    if ( ufun == 1 )
      dxdt = @(t,x) f1(t,x,@ufun1);
    else
      dxdt = @(t,x) f1(t,x,@ufun2);
    end
    tspan = [ t0, t1 ];
    x0 = [ 1.0, 2.0 ];
    [ t, x ] = ode45 ( dxdt, tspan, x0 );

    plot ( t, x, 'linewidth', 3 );
    grid ( 'on' );
    xlabel ( '<-- t -->' );
    ylabel ( '<-- X(t) -->' );
    s = sprintf ( 'State variable X(t) for control #%d', ufun );
    title ( s );
    legend ( '{x_1}', '{x_2}' );
    filename = sprintf ( 'x_u%d.png', ufun );
    print ( '-dpng', filename );
    fprintf ( 1, '  Graphics saved as "%s"\n', filename );
%
%  Evaluate cost.
%
    cost = trapz ( t, log(x(:,2)) );
    fprintf ( '  Cost functional = %g\n', cost );
%
%  Integrate backward in time for LAMBDA.
%
    t0 = 0.0;
    t1 = 5.0;
    if ( ufun == 1 )
      dlamdt = @(t,lam) f2(t,lam,x,@ufun1);
    else
      dlamdt = @(t,lam) f2(t,lam,x,@ufun2);
    end
    tspan = [ t1, t0 ];
    lam0 = [ 0.0, 0.0 ];
    [ t, lam ] = ode45 ( dlamdt, tspan, lam0 );

    plot ( t, lam, 'linewidth', 3 );
    grid ( 'on' );
    xlabel ( '<-- t -->' );
    ylabel ( '<-- {\lambda}(t) -->' );
    s = sprintf ( 'Adjoint variable {\\lambda}(t) for control #%d', ufun );
    title ( s );
    legend ( '{\lambda_1}', '{\lambda_2}' );
    filename = sprintf ( 'lam_u%d.png', ufun );
    print ( '-dpng', filename );
    fprintf ( 1, '  Graphics saved as "%s"\n', filename );

   end

  return
end
function dxdt = f1 ( t, x, ufun )
%
%  Discussion:
%
%    The text has
%      x1' = u x1
%      x2' = ( 1-u) x2
%    but apparently the correct form is
%      x1' = u x1
%      x2' = ( 1-u) x1
%
  dxdt = [         ufun(t)   * x(1); ...
           ( 1.0 - ufun(t) ) * x(1) ];    
  return
end
function value = ufun1 ( t )

  value = 0.0 * ( 0.0 <= t & t <  1.0 ) ...
        + 0.5 * ( 1.0 <= t & t <  3.0 ) ...
        + 1.0 * ( 3.0 <= t & t <= 5.0 );

  return
end
function value = ufun2 ( t )

  value = 1.0 * ( 0.0 <= t & t <  2.0 ) ...
        + 0.5 * ( 2.0 <= t & t <  4.0 ) ...
        + 0.0 * ( 4.0 <= t & t <= 5.0 );

  return
end
function dlamdt = f2 ( t, lam, x, ufun )

  dlamdt = [ - lam(1) * ufun(t) - lam(2) * ( 1.0 - ufun(t) ); ...
             - 1.0 / x(2)  ];
  return
end
