/*                                                                            */
/* 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                                                         */
/*                                                                            */
/* Portions Copyright (c) Year, Organization                                  */
/* All rights reserved.                                                       */
/*                                                                            */
/* Contributors:                                                              */
/*    Anshul Chawla                                                           */
/*     
#include <math.h>

/*                                                                            */
/* Auxiliary files for the KDS 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 beta = params[PARAM_beta];
  double const gamma = params[PARAM_gamma];
  double const A = params[PARAM_A];
  double const R = params[PARAM_R];

  double fc;

  if (Rij <= R)
  {
    fc = A;
  }
  else
  {
    fc = A*exp(-beta * pow(Rij - R, gamma));
  }

  return fc;
}

static double dfc_dR(double const * const params, double const Rij)
{
  /* Unpack parameters */

  double const beta = params[PARAM_beta];
  double const gamma = params[PARAM_gamma];
  double const A = params[PARAM_A];
  double const R = params[PARAM_R];

  double dfc_dR;
  dfc_dR = 0.0;

  if (Rij <= R)
  {
      dfc_dR = 0.0;
  }
  else
  {
      dfc_dR = -A*beta * gamma * pow(Rij - R, gamma - 1.0) * exp(-beta * pow(Rij - R, gamma));
  }

  return dfc_dR;
}

static double fR(double const * params, double const Rij)
{
  /* Unpack parameters */

  double const theta = params[PARAM_theta];

  return exp(-theta * Rij);

}

static double dfR_dRij(double const * const params, double const Rij)
{
  /* Unpack parameters */
  double const theta = params[PARAM_theta];

  return  -theta * exp(-theta * Rij);
}

static double fA(double const * const params, double const Rij)
{
  /* Unpack parameters */
  double const lambda = params[PARAM_lambda];

  return  -exp(-lambda * Rij);
}

static double dfA_dRij(double const * const params, double const Rij)
{
  /* Unpack parameters */
  double const lambda = params[PARAM_lambda];

  return lambda * exp(-lambda * Rij);
}

/*                                                                            */
/* Define functions used in three-body calculations                           */
/*                                                                            */
static double g(double const * const params, double const theta_jik)
{
  /* Unpack parameters */

  double const eta = params[PARAM_eta];
  double const thetanot = params[PARAM_thetanot];
  double const diff = theta_jik - thetanot;
  
  return (cos(eta*diff) - 1.0);
}

static double dg_dcostheta(double const * const params, double const theta_jik)
{
  /* Unpack parameters */

  double const eta = params[PARAM_eta];
  double const thetanot = params[PARAM_thetanot];
  double const diff = theta_jik - thetanot;
  
  return (((sin(eta*diff) * eta )/(sin(theta_jik))));
}