function tub ( ps, N )

%*****************************************************************************80
%
%% tub() optimizes the allocation of search boats for the Unsinkable Tub.
%
%  Discussion:
%
%    This code computes the optimal allocation of search boats
%    looking for the UNSINKABLE TUB.
%    The tub is sinking either 10 miles south (probability P1) or 
%    10 miles north of the station (probability P2 = 1 - P1).
%    N search boats are available, and each has a PS probability of locating
%    the tub.  How can we maximize the chance of locating the tub?
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    28 January 2018
%
%  Author:
%
%    Paul Nahin
%
%  Reference:
%
%    Paul Nahin,
%    Dueling Idiots and Other Probability Puzzlers,
%    Princeton, 2012,
%    ISBN: 978-0691155005.
%
%  Input:
%
%    real PS, the probability a search boat will find the tub.
%
%    integer N, the number of search boats available.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'tub():\n' );
  fprintf ( 1, '  Compute the optimal allocation of search boats\n' );
  fprintf ( 1, '  looking for the UNSINKABLE TUB.\n' );
  fprintf ( 1, '  The tub is sinking either 10 miles south (probability P1) or \n' );
  fprintf ( 1, '  10 miles north of the station (probability P2 = 1 - P1).\n' );
  fprintf ( 1, '  N search boats are available, and each has a PS probability of locating\n' );
  fprintf ( 1, '  the tub.  How can we maximize the chance of locating the tub?\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Report results graphically.\n' );

  for index = 1 : 999

    p1 = index / 1000;
    p2 = 1 - p1;
    n = ( N + ( log ( p2 / p1 ) / log ( 1 - ps ) ) ) / 2;
    nl = floor ( n );
    nu = nl + 1;
    Pl = p1 * ( 1 - ( 1 - ps ) ^ nl ) + p2 * ( 1 - ( 1 - ps ) ^ ( N - nl ) );
    Pu = p1 * ( 1 - ( 1 - ps ) ^ nu ) + p2 * ( 1 - ( 1 - ps ) ^ ( N - nu ) );

    n = min ( n, N );
    n = max ( n, 0 );

    if ( 0 < n && n < N )
      if ( Pu < Pl )
        n = nl;
      else
        n = nu;
      end
    end

    answer(index) = n;

  end
%
%  Graphics.
%
  clf ( );
  p = 0.001 : 0.001 : 0.999;
  plot ( p, answer, 'linewidth', 3 );
  grid ( 'on' );
  ylabel ( 'n, Search boats assigned to Island #1');
  xlabel ( 'p1, probability the UNSINKABLE TUB is at Island #1' );
  title ( 'The Unsinkable Tub' );
  filename = 'tub.png';
  print ( '-dpng', filename );
  fprintf ( '  Graphics saved as "%s"\n', filename );

  return
end
