/*                                                                            */
/* CDDL HEADER START                                                          */
/*                                                                            */
/* The contents of this file are subject to the terms of the Common           */
/* Development and Distribution License Version 1.0 (the "License").          */
/*                                                                            */
/* You can obtain a copy of the license at                                    */
/* http://www.opensource.org/licenses/CDDL-1.0.  See the License for the      */
/* specific language governing permissions and limitations under the License. */
/*                                                                            */
/* When distributing Covered Code, include this CDDL HEADER in each file and  */
/* include the License file in a prominent location with the name             */
/* LICENSE.CDDL.  If applicable, add the following below this CDDL HEADER,    */
/* with the fields enclosed by brackets "[]" replaced with your own           */
/* identifying information:                                                   */
/*                                                                            */
/* Portions Copyright (c) [yyyy] [name of copyright owner].                   */
/* All rights reserved.                                                       */
/*                                                                            */
/* CDDL HEADER END                                                            */
/*                                                                            */

/*                                                                            */
/* Portions Copyright (c) 2019, Regents of the University of Minnesota.       */
/* All rights reserved.                                                       */
/*                                                                            */
/* Contributors:                                                              */
/*    Daniel S. Karls                                                         */
/*    Ellad B. Tadmor                                                         */
/*    Ryan S. Elliott                                                         */
/*                                                                            */
/* <REPLACE: Add copyright information as needed.>                            */
/* Portions Copyright (c) Year, Organization                                  */
/* All rights reserved.                                                       */
/* </REPLACE>                                                                 */
/*                                                                            */
/* <REPLACE: Add contributors as needed.>                                     */
/* Contributors:                                                              */
/* </REPLACE>                                                                 */
/*                                                                            */

#include <math.h>

/*                                                                            */
/* Auxiliary files for the Tersoff T2 model driver                            */
/*                                                                            */
/* Functions in this file are only used by bondorder.inc. They are not        */
/* required by ThreeBodyBondOrder.c.                                          */
/*                                                                            */

static double const PI = 3.141592653589793;

/*                                                                            */
/* Define functions used in two-body calculations */
/*                                                                            */
static double fc(double const * const params, double const Rij)
{
  /* Unpack parameters */
  double const R = params[PARAM_R];
  double const D = params[PARAM_D];

  double fc;

  if (Rij <= R - D) { fc = 1.0; }
  else if (R - D < Rij && Rij < R + D)
  {
    fc = 0.5 * (1 - sin(0.5 * PI * (Rij - R) / D));
  }
  else
  {
    fc = 0.0;
  }

  return fc;
}

static double dfc_dR(double const * const params, double const Rij)
{
  /* Unpack parameters */
  double const R = params[PARAM_R];
  double const D = params[PARAM_D];

  double dfc_dR;

  dfc_dR = 0.0;

  if (Rij > (R - D) && Rij < (R + D))
  { dfc_dR = -(0.25 * PI / D) * cos(0.5 * PI * (Rij - R) / D); }

  return dfc_dR;
}

static double fR(double const * params, double const Rij)
{
  /* Unpack parameters */
  double const A = params[PARAM_A];
  double const lambda1 = params[PARAM_lambda1];

  return A * exp(-lambda1 * Rij);
}

static double dfR_dRij(double const * const params, double const Rij)
{
  /* Unpack parameters */
  double const A = params[PARAM_A];
  double const lambda1 = params[PARAM_lambda1];

  return -A * lambda1 * exp(-lambda1 * Rij);
}

static double fA(double const * const params, double const Rij)
{
  /* Unpack parameters */
  double const B = params[PARAM_B];
  double const lambda2 = params[PARAM_lambda2];

  return -B * exp(-lambda2 * Rij);
}

static double dfA_dRij(double const * const params, double const Rij)
{
  /* Unpack parameters */
  double const B = params[PARAM_B];
  double const lambda2 = params[PARAM_lambda2];

  return B * lambda2 * exp(-lambda2 * Rij);
}

/*                                                                            */
/* Define functions used in three-body calculations                           */
/*                                                                            */
static double g(double const * const params, double const costheta_jik)
{
  /* Unpack parameters */
  double const c = params[PARAM_c];
  double const d = params[PARAM_d];
  double const h = params[PARAM_h];

  return 1.0 + (c * c) / (d * d)
         - (c * c) / (d * d + (h - costheta_jik) * (h - costheta_jik));
}

static double dg_dcostheta(double const * const params, double const costheta_jik)
{
  /* Unpack parameters */
  double const c = params[PARAM_c];
  double const d = params[PARAM_d];
  double const h = params[PARAM_h];

  double denom;

  denom = ((d * d) + (h - costheta_jik) * (h - costheta_jik));
  denom *= denom;

  return -2 * (c * c) * (h - costheta_jik) / denom;
}