# include <stdlib.h>
# include <stdio.h>
# include <time.h>

# include "piecewise_linear_product_integral.h"

/******************************************************************************/

double piecewise_linear_product_integral ( double a, double b, int f_num, 
  double f_x[], double f_v[], int g_num, double g_x[], double g_v[] )

/******************************************************************************/
/*
  Purpose:

    PIECEWISE_LINEAR_PRODUCT_INTEGRAL: piecewise linear product integral.

  Discussion:

    We are given two piecewise linear functions F(X) and G(X) and we wish
    to compute the exact value of the integral

      INTEGRAL = Integral ( A <= X <= B ) F(X) * G(X) dx

    The functions F(X) and G(X) are defined as tables of coordinates X and
    values V.  A piecewise linear function is evaluated at a point X by 
    evaluating the interpolant to the data at the endpoints of the interval 
    containing X.  

    It must be the case that A <= B.

    It must be the case that the node coordinates F_X(*) and G_X(*) are
    given in ascending order.

    It must be the case that:

      F_X(1) <= A and B <= F_X(F_NUM)
      G_X(1) <= A and B <= G_X(G_NUM)

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    04 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, double A, B, the limits of integration.

    Input, int F_NUM, the number of nodes for F.

    Input, double F_X[F_NUM], the node coordinates for F.

    Input, double F_V[F_NUM], the nodal values for F.

    Input, int G_NUM, the number of nodes for G.

    Input, double G_X[G_NUM], the node coordinates for G.

    Input, double G_V[G_NUM], the nodal values for G.

    Output, double INTEGRAL, the integral of F(X) * G(X)
    from A to B.
*/
{
  double bit;
  int f_left;
  double f0;
  double f1;
  double fl;
  double fr;
  int g_left;
  double g0;
  double g1;
  double gl;
  double gr;
  double h0;
  double h1;
  double h2;
  int i;
  double integral;
  double xl;
  double xr;
  double xr_max;

  integral = 0.0;

  if ( f_x[f_num-1] <= a || g_x[g_num-1] <= a )
  {
    return integral;
  }

  if ( f_num < 2 || g_num < 2 )
  {
    return integral;
  }

  xr = a;

  f_left = 0;
  r8vec_bracket3 ( f_num, f_x, xr, &f_left );
  fr = f_v[f_left] + ( xr - f_x[f_left] ) * ( f_v[f_left+1] - f_v[f_left] ) 
    / ( f_x[f_left+1] - f_x[f_left] );

  g_left = 0;
  r8vec_bracket3 ( g_num, g_x, xr, &g_left );
  gr = g_v[g_left] + ( xr - g_x[g_left] ) * ( g_v[g_left+1] - g_v[g_left] ) 
    / ( g_x[g_left+1] - g_x[g_left] );

  xr_max = b;
  xr_max = r8_min ( xr_max, f_x[f_num-1] );
  xr_max = r8_min ( xr_max, g_x[g_num-1] );

  while ( xr < xr_max )
  {
/*
  Shift right values to left.
*/
    xl = xr;
    fl = fr;
    gl = gr;
/*
  Determine the new right values.
  The hard part is figuring out how to advance XR some, but not too much.
*/
    xr = xr_max;

    for ( i = 1; i <= 2; i++ )
    {
      if ( f_left + i <= f_num - 1 )
      {
        if ( xl < f_x[f_left+i] && f_x[f_left+i] < xr )
        {
          xr = f_x[f_left+i];
          break;
        }
      }
    }

    for ( i = 1; i <= 2; i++ )
    {
      if ( g_left + i <= g_num - 1 )
      {
        if ( xl < g_x[g_left+i] && g_x[g_left+i] < xr )
        {
          xr = g_x[g_left+i];
          break;
        }
      }
    }

    r8vec_bracket3 ( f_num, f_x, xr, &f_left );
    fr = f_v[f_left] + ( xr - f_x[f_left] ) * ( f_v[f_left+1] - f_v[f_left] ) 
      / ( f_x[f_left+1] - f_x[f_left] );

    r8vec_bracket3 ( g_num, g_x, xr, &g_left );
    gr = g_v[g_left] + ( xr - g_x[g_left] ) * ( g_v[g_left+1] - g_v[g_left] ) 
      / ( g_x[g_left+1] - g_x[g_left] );
/*
  Form the linear polynomials for F(X) and G(X) over [XL,XR],
  then the product H(X), integrate H(X) and add to the running total.
*/
    if ( r8_epsilon ( ) <= r8_abs ( xr - xl ) )
    {
      f1 = fl - fr;
      f0 = fr * xl - fl * xr;

      g1 = gl - gr;
      g0 = gr * xl - gl * xr;

      h2 = f1 * g1;
      h1 = f1 * g0 + f0 * g1;
      h0 = f0 * g0;

      h2 = h2 / 3.0;
      h1 = h1 / 2.0;

      bit = ( ( h2 * xr + h1 ) * xr + h0 ) * xr 
          - ( ( h2 * xl + h1 ) * xl + h0 ) * xl;

      integral = integral + bit / ( xr - xl ) / ( xr - xl );
    }
  }

  return integral;
}
/******************************************************************************/

double piecewise_linear_product_quad ( double a, double b, int f_num, 
  double f_x[], double f_v[], int g_num, double g_x[], double g_v[], 
  int quad_num )

/******************************************************************************/
/*
  Purpose:

    PIECEWISE_LINEAR_PRODUCT_QUAD: estimate piecewise linear product integral.

  Discussion:

    We are given two piecewise linear functions F(X) and G(X) and we wish
    to estimate the value of the integral

      INTEGRAL = Integral ( A <= X <= B ) F(X) * G(X) dx

    The functions F(X) and G(X) are defined as tables of coordinates X and
    values V.  A piecewise linear function is evaluated at a point X by 
    evaluating the interpolant to the data at the endpoints of the interval 
    containing X.  

    It must be the case that A <= B.

    It must be the case that the node coordinates F_X(*) and G_X(*) are
    given in ascending order.

    It must be the case that:

      F_X(1) <= A and B <= F_X(F_NUM)
      G_X(1) <= A and B <= G_X(G_NUM)

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    04 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, double A, B, the limits of integration.

    Input, int F_NUM, the number of nodes for F.

    Input, double F_X[F_NUM], the node coordinates for F.

    Input, double F_V[F_NUM], the nodal values for F.

    Input, int G_NUM, the number of nodes for G.

    Input, double G_X[G_NUM], the node coordinates for G.

    Input, double G_V[G_NUM], the nodal values for G.

    Input, int QUAD_NUM, the number of quadrature points.

    Output, double PIECEWISE_LINEAR_PRODUCT_QUAD, an estimate for the integral 
    of F(X) * G(X) from A to B.
*/
{
  double a2;
  double b2;
  int f_left;
  double fq;
  int g_left;
  double gq;
  int i;
  double quad;
  double xq;

  quad = 0.0;

  f_left = 0;
  g_left = 0;

  a2 = a;
  a2 = r8_max ( a2, f_x[0] );
  a2 = r8_max ( a2, g_x[0] );

  b2 = b;
  b2 = r8_min ( b2, f_x[f_num-1] );
  b2 = r8_min ( b2, g_x[g_num-1] );

  for ( i = 1; i <= quad_num; i++ )
  {
    xq =  ( ( double ) (                2 * i - 1 ) * b2 
          + ( double ) ( 2 * quad_num - 2 * i + 1 ) * a2 )  
          / ( double ) ( 2 * quad_num             );

    r8vec_bracket3 ( f_num, f_x, xq, &f_left );

    fq = f_v[f_left] + ( xq - f_x[f_left] ) * ( f_v[f_left+1] - f_v[f_left] )
      / ( f_x[f_left+1] - f_x[f_left] );

    r8vec_bracket3 ( g_num, g_x, xq, &g_left );

    gq = g_v[g_left] + ( xq - g_x[g_left] ) * ( g_v[g_left+1] - g_v[g_left] ) 
      / ( g_x[g_left+1] - g_x[g_left] );

    quad = quad + fq * gq;
  }

  quad = quad * ( b - a ) / ( double ) ( quad_num );

  return quad;
}
/******************************************************************************/

double r8_abs ( double x )

/******************************************************************************/
/*
  Purpose:

    R8_ABS returns the absolute value of an R8.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    07 May 2006

  Author:

    John Burkardt

  Parameters:

    Input, double X, the quantity whose absolute value is desired.

    Output, double R8_ABS, the absolute value of X.
*/
{
  double value;

  if ( 0.0 <= x )
  {
    value = + x;
  }
  else
  {
    value = - x;
  }
  return value;
}
/******************************************************************************/

double r8_epsilon ( )

/******************************************************************************/
/*
  Purpose:

    R8_EPSILON returns the R8 round off unit.

  Discussion:

    R8_EPSILON is a number R which is a power of 2 with the property that,
    to the precision of the computer's arithmetic,
      1 < 1 + R
    but
      1 = ( 1 + R / 2 )

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    01 September 2012

  Author:

    John Burkardt

  Parameters:

    Output, double R8_EPSILON, the R8 round-off unit.
*/
{
  const double value = 2.220446049250313E-016;

  return value;
}
/******************************************************************************/

double r8_max ( double x, double y )

/******************************************************************************/
/*
  Purpose:

    R8_MAX returns the maximum of two R8's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    07 May 2006

  Author:

    John Burkardt

  Parameters:

    Input, double X, Y, the quantities to compare.

    Output, double R8_MAX, the maximum of X and Y.
*/
{
  double value;

  if ( y < x )
  {
    value = x;
  }
  else
  {
    value = y;
  }
  return value;
}
/******************************************************************************/

double r8_min ( double x, double y )

/******************************************************************************/
/*
  Purpose:

    R8_MIN returns the minimum of two R8's.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    07 May 2006

  Author:

    John Burkardt

  Parameters:

    Input, double X, Y, the quantities to compare.

    Output, double R8_MIN, the minimum of X and Y.
*/
{
  double value;

  if ( y < x )
  {
    value = y;
  }
  else
  {
    value = x;
  }
  return value;
}
/******************************************************************************/

void r8vec_bracket3 ( int n, double t[], double tval, int *left )

/******************************************************************************/
/*
  Purpose:

    R8VEC_BRACKET3 finds the interval containing or nearest a given value.

  Discussion:

    An R8VEC is a vector of R8's.

    The routine always returns the index LEFT of the sorted array
    T with the property that either
    *  T is contained in the interval [ T[LEFT], T[LEFT+1] ], or
    *  T < T[LEFT] = T[0], or
    *  T > T[LEFT+1] = T[N-1].

    The routine is useful for interpolation problems, where
    the abscissa must be located within an interval of data
    abscissas for interpolation, or the "nearest" interval
    to the (extreme) abscissa must be found so that extrapolation
    can be carried out.

    This version of the function has been revised so that the value of
    LEFT that is returned uses the 0-based indexing natural to C++.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    31 May 2009

  Author:

    John Burkardt

  Parameters:

    Input, int N, length of the input array.

    Input, double T[N], an array that has been sorted into ascending order.

    Input, double TVAL, a value to be bracketed by entries of T.

    Input/output, int *LEFT.
    On input, if 0 <= LEFT <= N-2, LEFT is taken as a suggestion for the
    interval [ T[LEFT-1] T[LEFT] ] in which TVAL lies.  This interval
    is searched first, followed by the appropriate interval to the left
    or right.  After that, a binary search is used.
    On output, LEFT is set so that the interval [ T[LEFT], T[LEFT+1] ]
    is the closest to TVAL; it either contains TVAL, or else TVAL
    lies outside the interval [ T[0], T[N-1] ].
*/
{
  int high;
  int low;
  int mid;
/*
  Check the input data.
*/
  if ( n < 2 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8VEC_BRACKET3 - Fatal error\n" );
    fprintf ( stderr, "  N must be at least 2.\n" );
    exit ( 1 );
  }
/*
  If *LEFT is not between 0 and N-2, set it to the middle value.
*/
  if ( *left < 0 || n - 2 < *left )
  {
    *left = ( n - 1 ) / 2;
  }
/*
  CASE 1: TVAL < T[*LEFT]:
  Search for TVAL in (T[I],T[I+1]), for I = 0 to *LEFT-1.
*/
  if ( tval < t[*left] )
  {
    if ( *left == 0 )
    {
      return;
    }
    else if ( *left == 1 )
    {
      *left = 0;
      return;
    }
    else if ( t[*left-1] <= tval )
    {
      *left = *left - 1;
      return;
    }
    else if ( tval <= t[1] )
    {
      *left = 0;
      return;
    }
/*
  ...Binary search for TVAL in (T[I],T[I+1]), for I = 1 to *LEFT-2.
*/
    low = 1;
    high = *left - 2;

    for ( ; ; )
    {
      if ( low == high )
      {
        *left = low;
        return;
      }

      mid = ( low + high + 1 ) / 2;

      if ( t[mid] <= tval )
      {
        low = mid;
      }
      else
      {
        high = mid - 1;
      }
    }
  }
/*
  CASE 2: T[*LEFT+1] < TVAL:
  Search for TVAL in (T[I],T[I+1]) for intervals I = *LEFT+1 to N-2.
*/
  else if ( t[*left+1] < tval )
  {
    if ( *left == n - 2 )
    {
      return;
    }
    else if ( *left == n - 3 )
    {
      *left = *left + 1;
      return;
    }
    else if ( tval <= t[*left+2] )
    {
      *left = *left + 1;
      return;
    }
    else if ( t[n-2] <= tval )
    {
      *left = n - 2;
      return;
    }
/*
  ...Binary search for TVAL in (T[I],T[I+1]) for intervals I = *LEFT+2 to N-3.
*/
    low = *left + 2;
    high = n - 3;

    for ( ; ; )
    {

      if ( low == high )
      {
        *left = low;
        return;
      }

      mid = ( low + high + 1 ) / 2;

      if ( t[mid] <= tval )
      {
        low = mid;
      }
      else
      {
        high = mid - 1;
      }
    }
  }
/*
  CASE 3: T[*LEFT] <= TVAL <= T[*LEFT+1]:
  T is just where the user said it might be.
*/
  else
  {
  }

  return;
}
/******************************************************************************/

void timestamp ( void )

/******************************************************************************/
/*
  Purpose:

    TIMESTAMP prints the current YMDHMS date as a time stamp.

  Example:

    31 May 2001 09:45:54 AM

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    24 September 2003

  Author:

    John Burkardt

  Parameters:

    None
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}