/*                                                                            */
/* 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 R1 = params[PARAM_R1];
  double const R2 = params[PARAM_R2];


  double fc;

  if (Rij <= R1) { fc = 1.0; }
  else if (R1 < Rij && Rij < R2)
  {
    fc = 0.5 + 9/16 * cos(PI * (Rij - R1) / (R2 - R1)) - 1 / 16 * cos(3 * PI * (Rij - R1) / (R2 - R1))  ;  
  }
  else
  {
    fc = 0.0;
  }

  return fc;
}

static double dfc_dR(double const * const params, double const Rij)
{
  /* Unpack parameters */
  double const R1 = params[PARAM_R1];
  double const R2 = params[PARAM_R2];

  double dfc_dR;

  dfc_dR = 0.0;

  if (Rij > R1 && Rij < R2)
  { dfc_dR = - 9 / 16 * PI / (R2 - R1) * sin(PI * (Rij - R1) / (R2 - R1)) + 1 / 16 * 3 * PI / (R2 - R1) * sin(3 * PI * (Rij - R1) / (R2 - R1));}

  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 h = params[PARAM_h];


  double const c1 = params[PARAM_c1];
  double const c2 = params[PARAM_c2];
  double const c3 = params[PARAM_c3];
  double const c4 = params[PARAM_c4];
  double const c5 = params[PARAM_c5];

  double cos_term;
  cos_term = (h - costheta_jik);

  return c1 + (c2 * cos_term * cos_term) / (c3 + cos_term * cos_term)
               * (1.0 + c4 * exp(-c5 * cos_term * cos_term));

}

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

  double costerm;

  double const c2 = params[PARAM_c2];
  double const c3 = params[PARAM_c3];
  double const c4 = params[PARAM_c4];
  double const c5 = params[PARAM_c5];
  
  double term1;
  double term2;
  double term3;

  costerm = (h - costheta_jik);
  term1 = 2 * c2 * c4 * c5 * costerm * costerm * costerm * exp(-c5 * costerm * costerm) / (c3 + costerm * costerm);
  term2 = 2 * c2 * costerm *costerm * costerm * (1.0 + c4 * exp(-c5 * costerm * costerm)) / (c3 + costerm * costerm) / (c3 + costerm * costerm);
  term3 = -2 * c2 * costerm * (1.0 + c4 * exp(-c5 * costerm * costerm)) / (c3 + costerm * costerm);

  return (term1 + term2 + term3) ;

}