function tub ( ps, N )

%*****************************************************************************80
%
%% TUB optimizes the allocation of search boats for the Unsinkable Tub.
%
%  Discussion:
%
%    This m-file 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?
%
%  Modified:
%
%    28 January 2018
%
%  Author:
%
%    Paul Nahin
%
%  Reference:
%
%    Paul Nahin,
%    Dueling Idiots and Other Probability Puzzlers,
%    Princeton, 2012,
%    ISBN: 978-0691155005.
%
%  Parameters:
%
%    Input, real PS, the probability a search boat will find the tub.
%
%    Input, integer N, the number of search boats available.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TUB:\n' );
  fprintf ( 1, '  This m-file computes 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 ( Pl > Pu )
        n = nl;
      else
        n = nu;
      end
    end

    answer(index) = n;

  end
%
%  Graphics.
%
  clf
  p = .001 : .001 : .999;
  plot ( p, answer )
  grid
  ylabel('n, Search boats assigned to Island #1')
  xlabel('p1, probability the UNSINKABLE TUB is at Island #1')
  title('Fig.16.3')
  print ( '-dpng', 'tub.png' );

  return
end
