# include <cmath>
# include <cstdlib>

using namespace std;

# include "laplace_radial_exact.hpp"

//****************************************************************************80

void laplace_radial_2d_exact ( int n, double x[], double y[], double a, 
  double b, double u[], double ux[], double uy[], double uxx[], double uxy[], 
  double uyy[] )

//****************************************************************************80
//
//  Purpose:
//
//    laplace_radial_2d_exact() evaluates an exact radial solution of the Laplace equation.
//
//  Discussion:
//
//    The equation is:
//      Uxx + Uyy = 0
//
//    The radial solution in 2D, with parameters a and b, is:
//      R = sqrt ( X^2 + Y^2 )
//      U(X,Y) = a * log ( R ) + b
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    04 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer n: the number of points
//
//    real X, Y: the coordinates of the points.
//
//    real a, b: parameters.
//
//  Output:
//
//    real U, Ux, Uy, Uxx, Uxy, Uyy: the solution and first and second
//    partial derivatives evaluated at the points (X,Y).
//
{
  int i;
  double r;

  for ( i = 0; i < n; i++ )
  {
    r = sqrt ( pow ( x[i], 2 ) + pow ( y[i], 2 ) );

    u[i] = a * log ( r ) + b;
    ux[i] = a * x[i] / ( pow ( x[i], 2 ) + pow ( y[i], 2 ) );
    uy[i] = a * y[i] / ( pow ( x[i], 2 ) + pow ( y[i], 2 ) );
    uxx[i] = a * ( - 2 * pow ( x[i], 2 ) 
      / ( pow ( x[i], 2 ) + pow ( y[i], 2 ) ) + 1 ) 
      / ( pow ( x[i], 2 ) + pow ( y[i], 2 ) );
    uxy[i] = - 2 * a * x[i] * y[i] 
      / pow ( ( pow ( x[i], 2 ) + pow ( y[i], 2 ) ), 2 );
    uyy[i] = a * ( - 2 * pow ( y[i], 2 ) 
      / ( pow ( x[i], 2 ) + pow ( y[i], 2 ) ) + 1 ) 
      / ( pow ( x[i], 2 ) + pow ( y[i], 2 ) );
  }

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

void laplace_radial_3d_exact ( int n, double x[], double y[], double z[], 
  double a, double b, double u[], double ux[], double uy[], double uz[], 
  double uxx[], double uxy[], double uxz[], double uyy[], double uyz[], 
  double uzz[] )

//****************************************************************************80
//
//  Purpose:
//
//
//    laplace_radial_3d_exact() evaluates an exact radial solution of the Laplace equation.
//
//  Discussion:
//
//    The equation is:
//      Uxx + Uyy = 0
//
//    The radial solution in 3D, with parameters a and b, is:
//      R = sqrt ( X^2 + Y^2 + Z^2 )
//      U(X,Y,Z) = a * log ( R ) + b
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    04 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer n: the number of points.
//
//    real X, Y, Z: the coordinates of the points.
//
//    real a, b: shape parameters.
//
//  Output:
//
//    real U, Ux, Uy, Uz, Uxx, Uxy, Uxz, Uyy, Uyz, Uzz: 
//    the solution and first and second partial derivatives evaluated 
//    at the points (X,Y).
//
{
  int i;
  double r;

  for ( i = 0; i < n; i++ )
  {
    r = sqrt ( pow ( x[i], 2 ) + pow ( y[i], 2 ) + pow ( z[i], 2 ) );

    u[i] = - a / r + b;
    ux[i] = a * x[i] 
      / pow ( pow ( x[i], 2 ) + pow ( y[i], 2 ) + pow ( z[i], 2 ), 1.5 );
    uy[i] = a * y[i] 
      / pow ( pow ( x[i], 2 ) + pow ( y[i], 2 ) + pow ( z[i], 2 ), 1.5 );
    uz[i] = a * z[i] 
      / pow ( pow ( x[i], 2 ) + pow ( y[i], 2 ) + pow ( z[i], 2 ), 1.5 );
    uxx[i] = a * ( - 3 * pow ( x[i], 2 ) 
      / ( pow ( x[i], 2 ) + pow ( y[i], 2 ) + pow ( z[i], 2 ) ) + 1 ) 
      / pow ( pow ( x[i], 2 ) + pow ( y[i], 2 ) + pow ( z[i], 2 ), 1.5 );
    uxy[i] = - 3 * a * x[i] * y[i] 
      / pow ( pow ( x[i], 2 ) + pow ( y[i], 2 ) + pow ( z[i], 2 ), 2.5 );
    uxz[i] = - 3 * a * x[i] * z[i] 
      / pow ( pow ( x[i], 2 ) + pow ( y[i], 2 ) + pow ( z[i], 2 ), 2.5 );
    uyy[i] = a * ( - 3 * pow ( y[i], 2 ) 
      / ( pow ( x[i], 2 ) + pow ( y[i], 2 ) + pow ( z[i], 2 ) ) + 1 ) 
      / pow ( pow ( x[i], 2 ) + pow ( y[i], 2 ) + pow ( z[i], 2 ), 1.5 );
    uyz[i] = - 3 * a * y[i] * z[i] 
      / pow ( pow ( x[i], 2 ) + pow ( y[i], 2 ) + pow ( z[i], 2 ), 2.5 );
    uzz[i] = a * ( - 3 * pow ( z[i], 2 ) 
      / ( pow ( x[i], 2 ) + pow ( y[i], 2 ) + pow ( z[i], 2 ) ) + 1 ) 
      / pow ( pow ( x[i], 2 ) + pow ( y[i], 2 ) + pow ( z[i], 2 ), 1.5 );
  }

  return;
}