function pdf = triangle_right_distance_pdf ( v, d )

%*****************************************************************************80
%
%% triangle_right_distance_pdf: PDF for the right triangle distance problem.
%
%  Discussion:
%
%    This calculation is for a right triangle with vertices given in V.
%
%    The triangle will be mapped to
%
%       A
%       |\
%     b |  \ c
%       |    \
%       C-----B
%          a
%
%    The PDF is 0 if D is less than 0, or greater than c.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2019
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Uwe Baesel,
%    The distribution function of the distance between two random points
%    in a right-angled triangle,
%    arXiv:1208.6228v2 17 September 2012.
%
%  Parameters:
%
%    Input, real V(3,2), the vertices of the right triangle.
%
%    Input, real D, the distance.
%
%    Output, real PDF, the probability density function for that distance.
%

%
%  Force D to be a column vector!
%
  d = d(:);
%
%  1: Check that the triangle is close enough to a right triangle.
%
  tol = 100.0 * sqrt ( eps ) * max ( max ( v.^2 ) );

  err = triangle_right_error ( v );

  if ( tol < err )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'triangle_right_distance_pdf - Fatal error!\n' );
    fprintf ( 1, '  Input triangle rightness error = %g\n', err );
    return
  end
%
%  2: Compute the triangle side lengths a <= b <= c.
%
  [ a, b, c ] = triangle_sides ( v );

  alpha = atan2 ( a, b );
  perimeter = a + b + c;
  area = 0.5 * a * b;
  h = a * cos ( alpha );
%
%  3: Evaluate the PDF formula.
%
  pdf = ( 2.0 * d / area ) ...
    .* ( pi + ( jstar ( d, a, b, c, alpha, h ) - perimeter * d ) / area );

  return
end
function value = hstar ( indx, t, a, b, c, alpha, h )

%*****************************************************************************80
%
%% HSTAR evaluates the H-star functions #1 through #4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 March 2019
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Uwe Baesel,
%    The distribution function of the distance between two random points
%    in a right-angled triangle,
%    arXiv:1208.6228v2 17 September 2012.
%
%  Parameters:
%
%    Input, integer INDX, the index of the function.
%
%    Input, real T, the arguments.
%
%    Input, real A, B, the length of the two legs of the triangle.  A <= B.
%
%    Input, real C, the length of the hypotenuse.
%
%    Input, real ALPHA, the angle BAC, in radians.
%
%    Input, real H, the value of A * cos ( ALPHA ).
%
%    Output, real VALUE, the value of the function.
%
  if ( indx == 1 )
    value = theta ( 1, a, b, alpha ) * t.^2 / 8.0;
  elseif ( indx == 2 )
    value = theta ( 2, a, b, alpha ) * t.^2 / 8.0 ...
      + theta ( 3, a, b, alpha ) * lstar ( 1, t, h ) ...
      + c * lstar ( 2, t, h );
  elseif ( indx == 3 )
    value = a * t + theta ( 4, a, b, alpha ) * t.^2 / 2.0 ...
      + 0.5 * theta ( 3, a, b, alpha ) * lstar ( 1, t, h ) ...
      + 0.5 * c * lstar ( 2, t, h ) ...
      + 0.5 * ( b / a ) * lstar ( 1, t, a ) ...
      + 0.5 * b * lstar ( 2, t, a );
  elseif ( indx == 4 )
    value = theta ( 5, a, b, alpha ) * t ...
      + theta ( 6, a, b, alpha ) * t.^2 / 8.0 ...
      + 0.5 * ( b / a ) * lstar ( 1, t, a ) ...
      + 0.5 * b * lstar ( 2, t, a ) ...
      + 0.5 * ( a / b ) * lstar ( 1, t, b ) ...
      + 0.5 * a * lstar ( 2, t, b );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HSTAR - Fatal error!\n' );
    fprintf ( 1, '  Called with INDX = %d, but only values 1:4 are legal.\n', indx );
    error ( 'HSTAR - Fatal error!' );
  end

  return
end
function value = jstar ( d, a, b, c, alpha, h )

%*****************************************************************************80
%
%% JSTAR evaluates the J-star function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2019
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Uwe Baesel,
%    The distribution function of the distance between two random points
%    in a right-angled triangle,
%    arXiv:1208.6228v2 17 September 2012.
%
%  Parameters:
%
%    Input, real D(*), the distance between two points in the triangle.
%
%    Input, real A, B, the length of the two legs of the triangle.  A <= B.
%
%    Input, real C, the length of the hypotenuse.
%
%    Input, real ALPHA, the angle BAC, in radians.
%
%    Input, real H, the value of A * cos ( ALPHA ).
%
%    Output, real VALUE(*), the value of J*(D).
%
  n = length ( d );
  value = zeros ( n, 1 );

  i1 = 0.0 <= d & d < h;
  i2 = h   <= d & d < a;
  i3 = a   <= d & d < b;
  i4 = b   <= d & d <= c;

  value(i1) = jstars ( 1, 0.0, d(i1), a, b, c, alpha, h );
  value(i2) = jstars ( 1, 0.0, h,     a, b, c, alpha, h ) ...
            + jstars ( 2, h,   d(i2), a, b, c, alpha, h );
  value(i3) = jstars ( 1, 0.0, h,     a, b, c, alpha, h ) ...
            + jstars ( 2, h,   a,     a, b, c, alpha, h ) ...
            + jstars ( 3, a,   d(i3), a, b, c, alpha, h );
  value(i4) = jstars ( 1, 0.0, h,     a, b, c, alpha, h ) ...
            + jstars ( 2, h,   a,     a, b, c, alpha, h ) ...
            + jstars ( 3, a,   b,     a, b, c, alpha, h ) ...
            + jstars ( 4, b,   d(i4), a, b, c, alpha, h );

  return
end
function value = jstars ( indx, s, t, a, b, c, alpha, h )

%*****************************************************************************80
%
%% JSTARS evaluates the J-star functions #1 through #4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2019
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Uwe Baesel,
%    The distribution function of the distance between two random points
%    in a right-angled triangle,
%    arXiv:1208.6228v2 17 September 2012.
%
%  Parameters:
%
%    Input, integer INDX, the index of the function.
%
%    Input, real S, T, the arguments.
%
%    Input, real A, B, the length of the two legs of the triangle.  A <= B.
%
%    Input, real C, the length of the hypotenuse.
%
%    Input, real ALPHA, the angle BAC, in radians.
%
%    Input, real H, the value of A * cos ( ALPHA ).
%
%    Output, real VALUE, the value of the function.
%
  value = hstar ( indx, t, a, b, c, alpha, h ) ...
        - hstar ( indx, s, a, b, c, alpha, h );

  return
end
function value = lstar ( indx, t, a )

%*****************************************************************************80
%
%% LSTAR evaluates the L-star functions #1 through #2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2019
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Uwe Baesel,
%    The distribution function of the distance between two random points
%    in a right-angled triangle,
%    arXiv:1208.6228v2 17 September 2012.
%
%  Parameters:
%
%    Input, integer INDX, the index of the function.
%
%    Input, real T, A, the arguments.
%
%    Output, real VALUE, the value of the function.
%
  if ( indx == 1 )
    value = 0.5 * ( a * sqrt ( t.^2 - a.^2 ) + t.^2 .* asin ( a ./ t ) );
  elseif ( indx == 2 )
    value = sqrt ( t.^2 - a.^2 ) + a * asin ( a ./ t );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'LSTAR - Fatal error!\n' );
    fprintf ( 1, '  Called with INDX = %d, but only values 1:2 are legal.\n', indx );
    error ( 'LSTAR - Fatal error!' );
  end

  return
end
function value = theta ( indx, a, b, alpha )

%*****************************************************************************80
%
%% THETA evaluates the theta auxilliary functions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 February 2019
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Uwe Baesel,
%    The distribution function of the distance between two random points
%    in a right-angled triangle,
%    arXiv:1208.6228v2 17 September 2012.
%
%  Parameters:
%
%    Input, integer INDX, the index of the desired theta function.
%    1 <= INDX <= 6.
%
%    Input, real A, B, the length of the two legs of the triangle.  A <= B.
%
%    Input, real ALPHA, the angle BAC, in radians.
%
  if ( indx == 1 )
    value = ( a / b ) * ( 2.0 * alpha + pi ) + ( 2.0 * b / a ) * ( pi - alpha ) + 6.0;
  elseif ( indx == 2 )
    value = ( a / b ) * ( 2.0 * alpha - pi ) - ( 2.0 * b / a ) * alpha + 6.0;
  elseif ( indx == 3 )
    value = ( a / b ) + ( b / a );
  elseif ( indx == 4 )
    value = 1.0 - ( b / a ) * alpha;
  elseif ( indx == 5 )
    value = a + b;
  elseif ( indx == 6 )
    value = ( a / b ) * ( 2.0 * alpha - pi ) - ( 2.0 * b / a ) * alpha + 2.0;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'THETA - Fatal error!\n' );
    fprintf ( 1, '  Called with INDX = %d, but only values 1:6 are legal.\n', indx );
    error ( 'THETA - Fatal error!' );
  end

  return
end
