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 GNU LGPL license.
#
#  Modified:
#
#    08 February 2023
#
#  Author:
#
#    Original MATLAB version by Paul Nahin.
#    This version by John Burkardt.
#
#  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.
#
  print ( '' )
  print ( 'tub():' )
  print ( '  Compute the optimal allocation of search boats' )
  print ( '  looking for the UNSINKABLE TUB.' )
  print ( '  The tub is sinking either 10 miles south (probability P1) or ' )
  print ( '  10 miles north of the station (probability P2 = 1 - P1).' )
  print ( '  N search boats are available, and each has a PS probability of locating' )
  print ( '  the tub.  How can we maximize the chance of locating the tub?' )
  print ( '' )

  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 and n < N )
      if ( Pu < Pl )
        n = nl
      else
        n = nu
      end
    end

    answer(index) = n

  end
#
#  Graphics.
#
  plt.clf ( )
  p = .001 : .001 : .999
  plt.plot ( p, answer )
  plt.grid ( True )
  plt.ylabel ( 'n, Search boats assigned to Island #1' )
  plt.xlabel ( 'p1, probability the UNSINKABLE TUB is at Island #1' )
  plt.title ( 'The Unsinkable Tub' )
  filename = 'tub.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "' + filename + '"' )
  plt.close ( )

  return
end
