/*                                                                            */
/* 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:                                                              */
/*    Kevin M. Schmalbach                                                     */
/*                                                                            */

#include <math.h>

/*                                                                            */
/* Auxiliary files for the Gong model driver		                      */
/*                                                                            */
/* Functions in this file are only used by cluster.inc. They are not required */
/* by ThreeBodyCluster.c.                                                     */
/*                                                                            */

/*                                                                            */
/* Define functions used in two-body calculations                             */
/*                                                                            */
static void f2_df2(double const * const params,
                   double const r,
                   double * const f2,
                   double * const df2_dr)
{
  /* Unpack parameters */
  double const A = params[PARAM_A];
  double const B = params[PARAM_B];
  double const a = params[PARAM_a];
  double const p = params[PARAM_p];
  double const q = params[PARAM_q];

  double const r_pow_minus_p = pow(r, -p);
  double const r_pow_minus_q = pow(r, -q);
  double const r_minus_a = (r - a);
  double const r_minus_a_sq = r_minus_a * r_minus_a;
  double const exp_inv_r_minus_a = exp(1.0 / r_minus_a);

  if (r < a)
  {
    *f2 = A * (B * r_pow_minus_p - r_pow_minus_q) * exp_inv_r_minus_a;

    if (df2_dr != NULL)
    {
      *df2_dr = (*f2) * (-1. / r_minus_a_sq)
                + A * exp_inv_r_minus_a
                      * (-p * B * r_pow_minus_p + q * r_pow_minus_q) / r;
    }
  }
  else
  {
    *f2 = 0.0;
    if (df2_dr != NULL) { *df2_dr = 0.0; }
  }
}

/*                                                                            */
/* Define functions used in three-body calculations                           */
/*                                                                            */
static double g(double const costheta,
		double const * const params)
{
  double const c0 = params[PARAM_c0];
  double const c1 = params[PARAM_c1];
  return (costheta + 1.0 / 3.0) * (costheta + 1.0 / 3.0) * ((costheta+c0)*(costheta+c0)+c1);
}

static double dg_dcostheta(double const costheta,
			   double const * const params)
{
  double const c0 = params[PARAM_c0];
  double const c1 = params[PARAM_c1];
  return 2.0 * (costheta + 1.0 / 3.0)*((costheta+c0)*(costheta+c0)+c1)+(costheta+1.0/3.0)*(costheta+1.0/3.0)*(2*(costheta+c0));
}

static double h(double const * const params,
                double const r,
                double const s,
                double const costheta)
{
  /* Unpack parameters */
  double const lambda1 = params[PARAM_lambda1];
  double const gamma = params[PARAM_gamma];
  double const a = params[PARAM_a];

  if (r < a && s < a)
  {
    return lambda1 * exp(gamma / (r - a) + gamma / (s - a)) * g(costheta,params);
  }
  else
  {
    return 0.0;
  }
}

static void dh_drdsdt(double const * const params,
                      double const r,
                      double const s,
                      double const costheta,
                      double * const dh_dr,
                      double * const dh_ds,
                      double * const dh_dt)
{
  /* Unpack parameters */
  double const lambda1 = params[PARAM_lambda1];
  double const gamma = params[PARAM_gamma];
  double const a = params[PARAM_a];

  double t;
  double gval;
  double dgdcos;
  double expterm;

  double dcostheta_dr;
  double dcostheta_ds;
  double dcostheta_dt;

  if (r < a && s < a)
  {
    /* Compute third distance using law of cosines */
    t = sqrt(r * r + s * s - 2 * r * s * costheta);

    dcostheta_dr = (r * r - s * s + t * t) / (2 * r * r * s);
    dcostheta_ds = (s * s - r * r + t * t) / (2 * r * s * s);
    dcostheta_dt = -t / (r * s);

    gval = g(costheta,params);
    expterm = exp(gamma / (r - a) + gamma / (s - a));
    dgdcos = dg_dcostheta(costheta,params);

    *dh_dr = lambda1
             * (expterm * dgdcos * dcostheta_dr
                + gval * expterm * (-gamma / ((r - a) * (r - a))));
    *dh_ds = lambda1
             * (expterm * dgdcos * dcostheta_ds
                + gval * expterm * (-gamma / ((s - a) * (s - a))));
    *dh_dt = lambda1 * expterm * dgdcos * dcostheta_dt;
  }
  else
  {
    *dh_dr = 0.0;
    *dh_ds = 0.0;
    *dh_dt = 0.0;
  }
}

static void f3_df3(double const * const params,
                   double const Rij,
                   double const Rik,
                   double const Rjk,
                   double * const f3,
                   double * const df3_dRij,
                   double * const df3_dRik,
                   double * const df3_dRjk)
{
  double costheta_jik;

  double dh1_dRij;
  double dh1_dRik;
  double dh1_dRjk;

  /* Law of cosines to get angle from distances */
  costheta_jik = (Rij * Rij + Rik * Rik - Rjk * Rjk) / (2.0 * Rij * Rik);

  *f3 = h(params, Rij, Rik, costheta_jik);

  if (df3_dRij != NULL)
  {
    dh_drdsdt(params, Rij, Rik, costheta_jik, &dh1_dRij, &dh1_dRik, &dh1_dRjk);

    *df3_dRij = dh1_dRij;
    *df3_dRik = dh1_dRik;
    *df3_dRjk = dh1_dRjk;
  }

  return;
}