# include <math.h>
# include <stdlib.h>

# include "laplace_radial_exact.h"

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

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[] )

/******************************************************************************/
/*
  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;
}
/******************************************************************************/

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[] )

/******************************************************************************/
/*
  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;
}