lookup.c

/*----------------------------------------------------------------------
  SerialReax - Reax Force Field Simulator

  Copyright (2010) Purdue University
  Hasan Metin Aktulga, haktulga@cs.purdue.edu
  Joseph Fogarty, jcfogart@mail.usf.edu
  Sagar Pandit, pandit@usf.edu
  Ananth Y Grama, ayg@cs.purdue.edu

  This program is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of
  the License, or (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  See the GNU General Public License for more details:
  <http://www.gnu.org/licenses/>.
  ----------------------------------------------------------------------*/

#include "lookup.h"

#include "tool_box.h"
#include "two_body_interactions.h"


void Make_Lookup_Table( real xmin, real xmax, int n,
        lookup_function f, lookup_table* t )
{
    int i;

    t->xmin = xmin;
    t->xmax = xmax;
    t->n = n;
    t->dx = (xmax - xmin) / (n - 1);
    t->inv_dx = 1.0 / t->dx;
    t->a = (n - 1) / (xmax - xmin);
    t->y = (real*) smalloc( n * sizeof(real),
            "Make_Lookup_Table::t->y" );

    for (i = 0; i < n; i++)
    {
        t->y[i] = f(i * t->dx + t->xmin);
    }

    // fprintf(stdout,"dx = %lf\n",t->dx);
    // for(i=0; i < n; i++)
    //   fprintf( stdout,"%d %lf %lf %lf\n",
    //            i, i/t->a+t->xmin, t->y[i], exp(i/t->a+t->xmin) );
}


/* Fills solution into x. Warning: will modify c and d! */
static void Tridiagonal_Solve( const real *a, const real *b,
        real *c, real *d, real *x, unsigned int n)
{
    int i;
    real id;

    /* Modify the coefficients. */
    c[0] /= b[0]; /* Division by zero risk. */
    d[0] /= b[0]; /* Division by zero would imply a singular matrix. */
    for (i = 1; i < n; i++)
    {
        id = (b[i] - c[i - 1] * a[i]); /* Division by zero risk. */
        c[i] /= id;         /* Last value calculated is redundant. */
        d[i] = (d[i] - d[i - 1] * a[i]) / id;
    }