# include <cmath>
# include <cstdlib>

using namespace std;

# include "minimal_surface_exact.hpp"

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

void minimal_surface_catenoid_exact ( int n, double X[], double Y[], double a, 
  double U[], double Ux[], double Uy[], double Uxx[], double Uxy[], 
  double Uyy[] )

//****************************************************************************80
//
//  Purpose:
//
//    minimal_surface_catenoid_exact() evaluates catenoid minimal surface solution U(X,Y).
//
//  Discussion:
//
//    The minimal surface equation describes a function with zero mean curvature.
//
//    The equation is:
//      (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0
//
//    The helicoid solution is:
//      U(X,Y) = acosh ( a sqrt(X^2+Y^2) ) / a
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John D Cook,
//    Closed-form minimal surface solutions,
//    https://www.johndcook.com/blog/2025/03/31/minimal-surface-solutions/
//    Posted 31 March 2025.
//
//  Input:
//
//    int n: the number of points.
//
//    double X(n), Y(n): the coordinates of the points.
//
//    double a: a shape parameter.
//
//  Output:
//
//    double U(n), Ux(n), Uy(n), Uxx(n), Uxy(n), Uyy(n): 
//    the solution, first and second partial derivatives evaluated at (X,Y).
//
{
  int i;
  double n1;
  double n2;

  for ( i = 0; i < n; i++ )
  {
    n1 = X[i] * X[i] + Y[i] * Y[i];
    n2 = sqrt ( n1 );

    U[i] = acosh ( a * n2 ) / a;
    Ux[i] = X[i] / ( n2 * sqrt ( a * a * n1 - 1.0 ) );
    Uy[i] = Y[i] / ( n2 * sqrt ( a * a * n1 - 1.0 ) );
    Uxx[i] = - ( a * a * X[i] * X[i] / ( a * a * n1 - 1.0 ) 
      + X[i] * X[i] / n1 - 1.0 ) / ( n2 * sqrt ( a * a * n1 - 1.0 ) );
    Uxy[i] = - X[i] * Y[i] * ( a * a / ( a * a * n1 - 1.0 ) 
      + 1.0 / n1 ) / ( n2 * sqrt ( a * a * n1 - 1.0 ) );
    Uyy[i] = - ( a * a * Y[i] * Y[i] / ( a * a * n1 - 1.0 ) 
      + Y[i] * Y[i] / n1 - 1.0 ) / ( n2 * sqrt ( a * a * n1 - 1.0 ) );
  }

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

double *minimal_surface_catenoid_residual ( int n, double X[], double Y[], 
  double a )

//****************************************************************************80
//
//  Purpose:
//
//    minimal_surface_catenoid_residual(): minimal surface residual for catenoid solution.
//
//  Discussion:
//
//    The equation is:
//      (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0
//
//    The catenoid solution is:
//      U(X,Y) = acosh ( a * sqrt ( X^2 + Y^2 ) ) / a
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int n: the number of points.
//
//    double X(n), Y(n): the coordinates of the points.
//
//    double a: a shape parameter.
//
//  Output:
//
//    double R(n): the residual evaluated at the points (X,Y).
//
{
  int i;
  double *R;
  double *U;
  double *Ux;
  double *Uxx;
  double *Uxy;
  double *Uy;
  double *Uyy;

  U = new double[n];
  Ux = new double[n];
  Uy = new double[n];
  Uxx = new double[n];
  Uxy = new double[n];
  Uyy = new double[n];

  minimal_surface_catenoid_exact ( n, X, Y, a, U, Ux, Uy, Uxx, Uxy, Uyy );

  R = new double[n];

  for ( i = 0; i < n; i++ )
  {
    R[i] = ( 1.0 + pow ( Ux[i], 2 ) ) * Uyy[i] 
      - 2.0 * Ux[i] * Uxy[i] * Uy[i]
      + ( 1.0 + pow ( Uy[i], 2 ) ) * Uxx[i];
  }

  delete []  U;
  delete []  Ux;
  delete []  Uy;
  delete []  Uxx;
  delete []  Uxy;
  delete []  Uyy;

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

void minimal_surface_helicoid_exact ( int n, double X[], double Y[], double U[], 
  double Ux[], double Uy[], double Uxx[], double Uxy[], double Uyy[] )

//****************************************************************************80
//
//  Purpose:
//
//    minimal_surface_helicoid_exact() evaluates helicoid minimal surface solution U(X,Y).
//
//  Discussion:
//
//    The minimal surface equation describes a function with zero mean curvature.
//
//    The equation is:
//      (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0
//
//    The helicoid solution is:
//      U(X,Y) = arctan ( X / Y )
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John D Cook,
//    Closed-form minimal surface solutions,
//    https://www.johndcook.com/blog/2025/03/31/minimal-surface-solutions/
//    Posted 31 March 2025.
//
//  Input:
//
//    int n: the number of points.
//
//    double X(n), Y(n): the coordinates of the points.
//
//  Output:
//
//    double U(n), Ux(n), Uy(n), Uxx(n), Uxy(n), Uyy(n): 
//    the solution, first and second partial derivatives evaluated at (X,Y).
//
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    U[i] = atan ( X[i] / Y[i] );
    Ux[i] = 1.0 / ( Y[i] * ( pow ( X[i], 2 ) / pow ( Y[i], 2 ) + 1.0 ) );
    Uy[i] = - X[i] / ( pow ( Y[i], 2 ) * ( pow ( X[i], 2 ) 
      / pow ( Y[i], 2 ) + 1.0 ) );
    Uxx[i] = - 2.0 * X[i] / ( pow ( Y[i], 3 ) 
      * pow ( ( pow ( X[i], 2 ) / pow ( Y[i], 2 ) + 1.0 ), 2 ) );
    Uxy[i] = ( 2.0 * pow ( X[i], 2 ) / ( pow ( Y[i], 2 ) 
      * ( pow ( X[i], 2 ) / pow ( Y[i], 2 ) + 1.0 ) ) - 1.0 ) 
      / ( pow ( Y[i], 2 ) * ( pow ( X[i], 2 ) / pow ( Y[i], 2 ) + 1.0 ) );
    Uyy[i] = 2.0 * X[i] * ( - pow ( X[i], 2 ) / ( pow ( Y[i], 2 ) 
      * ( pow ( X[i], 2 ) / pow ( Y[i], 2 ) + 1.0 ) ) + 1.0 ) 
      / ( pow ( Y[i], 3 ) * ( pow ( X[i], 2 ) / pow ( Y[i], 2 ) + 1.0 ) );
  }

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

double *minimal_surface_helicoid_residual ( int n, double X[], double Y[] )

//****************************************************************************80
//
//  Purpose:
//
//    minimal_surface_helicoid_residual(): minimal surface residual for helicoid solution.
//
//  Discussion:
//
//    The equation is:
//      (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0
//
//    The helicoid solution is:
//      U(X,Y) = arctan ( X / Y )
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int n: the number of points.
//
//    double X(n), Y(n): the coordinates of the points.
//
//  Output:
//
//    double R(n): the residual evaluated at the points (X,Y).
//
{
  int i;
  double *R;
  double *U;
  double *Ux;
  double *Uxx;
  double *Uxy;
  double *Uy;
  double *Uyy;

  U = new double[n];
  Ux = new double[n];
  Uy = new double[n];
  Uxx = new double[n];
  Uxy = new double[n];
  Uyy = new double[n];

  minimal_surface_helicoid_exact ( n, X, Y, U, Ux, Uy, Uxx, Uxy, Uyy );

  R = new double[n];

  for ( i = 0; i < n; i++ )
  {
    R[i] = ( 1.0 + pow ( Ux[i], 2 ) ) * Uyy[i] 
      - 2.0 * Ux[i] * Uxy[i] * Uy[i] 
      + ( 1.0 + pow ( Uy[i], 2 ) ) * Uxx[i];
  }

  delete []  U;
  delete []  Ux;
  delete []  Uy;
  delete []  Uxx;
  delete []  Uxy;
  delete []  Uyy;

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

void minimal_surface_linear_exact ( int n, double X[], double Y[], double a, 
  double b, double c, double U[], double Ux[], double Uy[], double Uxx[], 
  double Uxy[], double Uyy[] )

//****************************************************************************80
//
//  Purpose:
//
//    minimal_surface_linear_exact() evaluates linear minimal surface solution U(X,Y).
//
//  Discussion:
//
//    The minimal surface equation describes a function with zero mean curvature.
//
//    The equation is:
//      (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0
//
//    The linear solution is:
//      U(X,Y) = a * X + b * Y + C
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John D Cook,
//    Closed-form minimal surface solutions,
//    https://www.johndcook.com/blog/2025/03/31/minimal-surface-solutions/
//    Posted 31 March 2025.
//
//  Input:
//
//    int n: the number of points.
//
//    double X(n), Y(n): the coordinates of the points.
//
//    double a, b, c: parameters.
//
//  Output:
//
//    double U(n), Ux(n), Uy(n), Uxx(n), Uxy(n), Uyy(n): 
//    the solution, first and second partial derivatives evaluated at (X,Y).
//
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    U[i] = a * X[i] + b * Y[i] + c;
    Ux[i] = a;
    Uy[i] = b;
    Uxx[i] = 0.0;
    Uxy[i] = 0.0;
    Uyy[i] = 0.0;
  }

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

double *minimal_surface_linear_residual ( int n, double X[], double Y[], 
  double a, double b, double c )

//****************************************************************************80
//
//  Purpose:
//
//    minimal_surface_linear_residual(): minimal surface residual for linear solution.
//
//  Discussion:
//
//    The equation is:
//      (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0
//
//    The linear solution is:
//      U(X,Y) = a * X + b * Y + C
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int n: the number of points.
//
//    double X(n), Y(n): the coordinates of the points.
//
//    double a, b, c: parameters.
//
//  Output:
//
//    double R(n): the residual evaluated at the points (X,Y).
//
{
  int i;
  double *R;
  double *U;
  double *Ux;
  double *Uxx;
  double *Uxy;
  double *Uy;
  double *Uyy;

  U = new double[n];
  Ux = new double[n];
  Uy = new double[n];
  Uxx = new double[n];
  Uxy = new double[n];
  Uyy = new double[n];

  minimal_surface_linear_exact ( n, X, Y, a, b, c, U, Ux, Uy, Uxx, Uxy, 
    Uyy );

  R = new double[n];

  for ( i = 0; i < n; i++ )
  {
    R[i] = ( 1.0 + pow ( Ux[i], 2 ) ) * Uyy[i] 
      - 2.0 * Ux[i] * Uxy[i] * Uy[i] 
      + ( 1.0 + pow ( Uy[i], 2 ) ) * Uxx[i];
  }

  delete []  U;
  delete []  Ux;
  delete []  Uy;
  delete []  Uxx;
  delete []  Uxy;
  delete []  Uyy;

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

void minimal_surface_scherk_exact ( int n, double X[], double Y[], double a, 
  double U[], double Ux[], double Uy[], double Uxx[], double Uxy[], 
  double Uyy[] )

//****************************************************************************80
//
//  Purpose:
//
//    minimal_surface_scherk_exact() evaluates Scherk minimal surface solution U(X,Y).
//
//  Discussion:
//
//    The minimal surface equation describes a function with zero mean curvature.
//
//    The equation is:
//      (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0
//
//    The Scherk solution is:
//      U(X,Y) = log ( cos ( a Y ) / cos ( a X ) ) / a
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    John D Cook,
//    Closed-form minimal surface solutions,
//    https://www.johndcook.com/blog/2025/03/31/minimal-surface-solutions/
//    Posted 31 March 2025.
//
//  Input:
//
//    int n: the number of points.
//
//    double X(n), Y(n): the coordinates of the points.
//
//    double a: a shape parameter.
//
//  Output:
//
//    double U(n), Ux(n), Uy(n), Uxx(n), Uxy(n), Uyy(n): 
//    the solution, first and second partial derivatives evaluated at (X,Y).
//
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    U[i] = log ( cos ( a * Y[i] ) / cos ( a * X[i] ) ) / a;
    Ux[i] =    sin ( a * X[i] ) / cos ( a * X[i] );
    Uy[i] =  - sin ( a * Y[i] ) / cos ( a * Y[i] );
    Uxx[i] =    a * ( pow ( sin ( a * X[i] ), 2 ) / 
                      pow ( cos ( a * X[i] ), 2 ) + 1.0 );
    Uxy[i] =  0.0;
    Uyy[i] =  - a * ( pow ( sin ( a * Y[i] ), 2 ) / 
                      pow ( cos ( a * Y[i] ), 2 ) + 1.0 );
  }

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

double *minimal_surface_scherk_residual ( int n, double X[], double Y[], 
  double a )

//****************************************************************************80
//
//  Purpose:
//
//    minimal_surface_scherk_residual(): minimal surface residual for Scherk solution.
//
//  Discussion:
//
//    The equation is:
//      (1+Ux^2) Uyy - 2 Ux Uy Uxy + (1+Uy^2) Uxx = 0
//
//    The Scherk solution is:
//      U(X,Y) = log ( cos ( a * Y ) / cos ( a * X ) ) / a
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 June 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int n: the number of points.
//
//    double X(n), Y(n): the coordinates of the points.
//
//    double a: a shape parameter.
//
//  Output:
//
//    double R(n): the residual evaluated at the points (X,Y).
//
{
  int i;
  double *R;
  double *U;
  double *Ux;
  double *Uxx;
  double *Uxy;
  double *Uy;
  double *Uyy;

  U = new double[n];
  Ux = new double[n];
  Uy = new double[n];
  Uxx = new double[n];
  Uxy = new double[n];
  Uyy = new double[n];

  minimal_surface_scherk_exact ( n, X, Y, a, U, Ux, Uy, Uxx, Uxy, Uyy );

  R = new double[n];

  for ( i = 0; i < n; i++ )
  {
    R[i] = ( 1.0 + pow ( Ux[i], 2 ) ) * Uyy[i]
      - 2.0 * Ux[i] * Uxy[i] * Uy[i] 
      + ( 1.0 + pow ( Uy[i], 2 ) ) * Uxx[i];
  }

  delete []  U;
  delete []  Ux;
  delete []  Uy;
  delete []  Uxx;
  delete []  Uxy;
  delete []  Uyy;

  return R;
}