xmd Code

/*
xmd - molecular dynamics for metals and ceramics
By Jonathan Rifkin <jon.rifkin@uconn.edu>
Copyright 1995-2004 Jonathan Rifkin

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.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/


/*
************************************************************************
History
************************************************************************
*/

/*
First pass at general function routine
*/


/*
************************************************************************
Compile Switches
************************************************************************
*/

/*  LOCAL:  Free standing compilation, does not use XMD routines  */
/*
#define LOCAL
*/




/*
************************************************************************
Include Files
************************************************************************
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "cdfunc.h"
#include "cdhouse.h"

#ifndef LOCAL
#include "iomngr.h"
#include "strsub.h"
#include "parse.h"
#endif



/*
************************************************************************
Defines
************************************************************************
*/
#define BOOLEAN int
#define TRUE 1
#define FALSE 0

#define X 0
#define Y 1
#define Z 2
#define NDIR 3

#define NBUFF 1024


#define DATA_NULL  0
#define DATA_TABLE 1
#define DATA_POLY  2


#define DERIV_ZERO 0
#define DERIV_ONE  1


#ifdef LOCAL
#define LIST FILE
#endif


/*
************************************************************************
Macros
************************************************************************
*/
#ifdef LOCAL

#define FGETS(STRING,SIZE,FILE) \
   fgets (STRING,SIZE,FILE)
#define DBLSTRF(STRING)  dblstr(STRING)
#define INTSTRF(STRING)  intstr(STRING)
#define STRCMPI(S1,S2)   strcasecmp(S1,S2)
#define STRNCMPI(S1,S2,N)  strncasecmp(S1,S2,N)

#else

#define FGETS(STRING,SIZE,FILE) \
   m_gets_list (STRING,SIZE,FILE)
#define DBLSTRF(STRING)  dblstrf(STRING)
#define INTSTRF(STRING)  intstrf(STRING)
#define STRCMPI(S1,S2)   strcmpi(S1,S2)
#define STRNCMPI(S1,S2)  strncmpi(S1,S2)

#endif

/*
************************************************************************
Type Definitions
************************************************************************
*/


/*
************************************************************************
Global Variables
************************************************************************
*/

#ifdef LOCAL
int Debug_g = TRUE;

FILE *inlist = stdin;
BOOLEAN CheckMem_g = FALSE;

#else

extern LIST         *inlist;

#endif

extern BOOLEAN UseCheckFunc_g;



/*
************************************************************************
Module-Wide Variables
************************************************************************
*/




/*
************************************************************************
Local Function Prototypes
************************************************************************
*/
static void ReadTableFunction      (Function_t *, char *, LIST *);
static void ReadPolyFunction       (Function_t *, char *, LIST *);
static void ScaleTableFunction     (Function_t *, double, int);
static void ScalePolyFunction      (Function_t *, double, int);
#if 0
static void CalcArrayTableFunction (Function_t *, char *, double *);
#endif

static void InterpolateTable       (Function_t *, double *);

static double GetTableValue (Function_t *, double);
static double GetTableDeriv (Function_t *, double);
static double GetNullValue (Function_t *F, double x);
static double GetPolyValue  (Function_t *, double);
static double GetPolyDeriv  (Function_t *, double);
static double OutOfRangeValue (Function_t *F, double x, int Deriv);

static double poly (double, int, double *);


#ifdef LOCAL
BOOLEAN IsWhiteSpace(char);
BOOLEAN IsBlank (char *);
char    *strhed (char **);
double  dblstr(char **);
double  intstr(char **);
void    CleanAfterError(void);
#endif



/*
************************************************************************
Exported Routines
************************************************************************
*/

int FUNC_Parse (Function_t *F, char *InputStr)
	{
	char *HeadStr = NULL;

	
	/*  Check commands other than input commands  */
	if (IsFirstToken("SCALE",InputStr))
		{
		strhed (&InputStr);
		FUNC_ParseScale (F, InputStr);
		return FUNC_COMMAND_SCALE;
		}

	/*  
	Pass on to FUNC_ParseRead
	*/
	else
		{
		FUNC_ParseRead (F, InputStr);
		return FUNC_COMMAND_READ;
		}
	}



void FUNC_ParseScale (Function_t *F, char *InputStr)
	{
	char *HeadStr = NULL;
	double Scale;

	HeadStr = strhed (&InputStr);
	Scale   = dblstrf(&InputStr);

	if (!STRCMPI("INPUT",HeadStr))
		{
		FUNC_Scale (F, Scale, FUNC_INPUT);
		}
	else if (!STRCMPI("OUTPUT",HeadStr))
		{
		FUNC_Scale (F, Scale, FUNC_OUTPUT);
		}
	else
		{
		ERROR_PREFIX
		printf ("Unrecognized option for SCALE command.\n");
		CleanAfterError();
		}
	}


/*
Read function table or polynomial from input
  - Input is throught macro FGETS(STRING,SIZE,FILE)
  - Received string of the form(s)
      [ TABLE ] n InteriorCutoff ExteriorCutoff
           x0 x1 ... xn-1
      POLY      n0 nn x0 ExteriorCutoff
           an0 an1 .. ann
      f(x) = an0*(x-xc)^n0 + ... ann*(x-xc)^nn
*/
void FUNC_ParseRead (Function_t *F, char *InputStr)
   {

   /*
   Test for valid function boundary
   */
   if (FUNC_EXTRAP_UNDEFINED==GET_INT_EXTRAP(F->BoundaryType))
      {
      printf ("INTERNAL ERROR: File %s Line %i\n",
         __FILE__, __LINE__);
      printf ("Invalid Interior Extrap value\n");
      CleanAfterError();
      }

   if (FUNC_EXTRAP_UNDEFINED==GET_EXT_EXTRAP(F->BoundaryType))
      {
      printf ("INTERNAL ERROR: File %s Line %i\n",
         __FILE__, __LINE__);
      printf ("Invalid Exterior Extrap value\n");
      CleanAfterError();
		}


   /*
   Polynomial Data Type
   */
	if (IsFirstToken("POLY",InputStr))
      {
		strhed (&InputStr);
      F->DataType = DATA_POLY;
      ReadPolyFunction (F, InputStr, inlist);
      }

   /*
   Table Data type (default)
   */
   else
      {
      F->DataType = DATA_TABLE;

      /*
      TABLE token found, remove it
      */
		if (IsFirstToken("TABLE",InputStr))
         {
			strhed (&InputStr);
         }

      ReadTableFunction (F, InputStr, inlist);
      }

   }


Function_t *FUNC_CreateNullFunction(void)
	{
	Function_t *F = NULL;

	/*  Allocate storage  */
	ALLOCATE (F, Function_t, 1)

	/*  Set function pointer to 'null' function  */
	F->f0 = (FunctionCall_t*) GetNullValue;
	F->f1 = (FunctionCall_t*) GetNullValue;
   F->DataType = DATA_NULL;

	return F;
	
	}

void FUNC_Free (Function_t **F)
	{
	if (*F!=NULL)
		{
		FREE ((*F)->Coeff)
		FREE ((*F))
		}
   }



/*  Scale function  */
void FUNC_Scale (Function_t *F, double Scale, int ScaleType)
   {

   /*  Check scale  */
      if (Scale==0 && FUNC_INPUT==ScaleType)
         {
         printf ("ERROR:  Cannot scale function input by 0.\n");
         CleanAfterError();
         }

   /*  Scale value common to both types  */
   if (FUNC_INPUT==ScaleType)
      {
      F->InteriorCutoff *= Scale;
      F->ExteriorCutoff *= Scale;
      F->InteriorSlope  /= Scale;
      F->ExteriorSlope  /= Scale;
      }
   else if (FUNC_OUTPUT==ScaleType)
      {
      F->InteriorValue *= Scale;
      F->ExteriorValue *= Scale;
      F->InteriorSlope *= Scale;
      F->ExteriorSlope *= Scale;
      }

   /*  Scale function specific values  */
   if (DATA_POLY==F->DataType)
      ScalePolyFunction(F, Scale, ScaleType);
   else
      ScaleTableFunction(F, Scale, ScaleType);
   }


double FUNC_GetValue (Function_t *F, double x)
   {
   if      (DATA_POLY ==F->DataType)
      return GetPolyValue (F, x);
	else if (DATA_TABLE==F->DataType)
      return GetTableValue (F, x);
	else
		return GetNullValue  (F, x);
   }



double FUNC_GetDeriv (Function_t *F, double x)
   {
   if (DATA_POLY==F->DataType)
      return GetPolyDeriv (F, x);
   else
      return GetTableDeriv (F, x);
   }




double FUNC_GetCutoff(Function_t *F)
   {
   return F->ExteriorCutoff;
   }



/*
************************************************************************
Local Functions
************************************************************************
*/
static void ReadTableFunction
(
Function_t *F,
char *InputStr,
LIST *InFile
)
	{
   int idata;
   int nread;
   int ErrorCode;
   char ReadBuffer[NBUFF];
   char TempBuffer[NBUFF];
	char *TokenStr = NULL;
	char *FormPtr  = NULL;
	char *TempPtr  = NULL;
	double *a = NULL;
	double *atmp = NULL;
	double DelR;
	double r;

	F->NumCoeff = INTSTRF(&InputStr) - 1;
	F->InteriorCutoff = DBLSTRF (&InputStr);
	F->ExteriorCutoff = DBLSTRF (&InputStr);
	F->f0 = (FunctionCall_t *) GetTableValue;
	F->f1 = (FunctionCall_t *) GetTableDeriv;

	/*  Sanity Check  */

   if (F->NumCoeff <  2)
      {
      printf
         ("ERROR:  ");
      printf
         ("There must be at least 2 data points in table potenial.\n");
		CleanAfterError();
      }

   if (F->ExteriorCutoff <= F->InteriorCutoff)
		{
      printf ( "ERROR:  Exterior Cutoff (%e)",
         F->ExteriorCutoff);
		printf (" is less than or equal to interior cutoff (%e).\n",
			F->InteriorCutoff);
      CleanAfterError();
      }

   F->TableFactor =
      F->NumCoeff / ( F->ExteriorCutoff - F->InteriorCutoff );

   /*  Make temporary space for input values  */
	ALLOCATE (atmp, double, F->NumCoeff+5);
   a     = atmp+2;

   /*  Make permanant space for coefficients  */
	ALLOCATE (F->Coeff, double, 6*F->NumCoeff)

	/*  Get next token  */
	TokenStr = strhed (&InputStr);

	/*  Check for valid token  */
	if (!IsBlank(TokenStr) && STRCMPI("FORMULA",TokenStr))
		{
		ERROR_PREFIX
		printf ("Unknown option \"%s\"\n", TokenStr);
		CleanAfterError();
		}

   /*  Read formula  */
	if (!STRCMPI("FORMULA",TokenStr))
      {

		FGETS(ReadBuffer,NBUFF,InFile);
		DelR = (F->ExteriorCutoff - F->InteriorCutoff)/F->NumCoeff;
		/*  
		Loop to NumCoeff+1 because, in addition to values at 
		start of all intervales, also include value at end of 
		of last interval.
		*/
      LOOP (idata, F->NumCoeff+1)
         {

         /*  Write R (in angstroms) to parsing variable R  */
         r = F->InteriorCutoff + idata * DelR;
         sprintf (TempBuffer, "R=%le", r);
         TempPtr = TempBuffer;
         dblstrf (&TempPtr);

         /*  Evaluate formula (will use value of R passed via dblstrf()  */
         FormPtr = ReadBuffer;
			a[idata] = evalform (&FormPtr, &ErrorCode);

         /*  Error checking  */
			if (ErrorCode)
            {
				printf
					(
					"ERROR (pr_read_potential):  Cannot parse formula\n"
					);
				printf
					(
					"   Parser returns following message: %s\n",
					parsemsg(ErrorCode)
					);
				CleanAfterError();
				}
			}
		}

	/*  Read input table  */
	else
		{
      nread = 0;
		while (nread<F->NumCoeff+1 && NULL!=FGETS(ReadBuffer,NBUFF,InFile))
			{
			TokenStr = ReadBuffer;
			while (*TokenStr && nread < F->NumCoeff+1)
				{
				a[nread] = DBLSTRF (&TokenStr);
				nread++;
            }
         }

      if (nread<F->NumCoeff+1)
         {
         printf ("ERROR: Potential table too short.\n");
			CleanAfterError();
         }
      }

   ASSERT (FUNC_EXTRAP_UNDEFINED!=GET_INT_EXTRAP(F->BoundaryType))
   ASSERT (FUNC_EXTRAP_UNDEFINED!=GET_EXT_EXTRAP(F->BoundaryType))

	/*  With data array in hand, call table interpolation routine  */
	InterpolateTable(F, a);

	/*  Free unneeded allocation  */
	FREE (atmp)
	}


#if 0
/*
Like ReadTableFunction() except data table passed as an array.
NOTE:  Array starts at a[-2] and runs to a[n+2]
Data stored in a[0] to a[n-1].
*/
static void CalcArrayTableFunction
(
Function_t *F,
char *InputStr,
double *a
)
	{

   F->NumCoeff = INTSTRF(&InputStr) - 1;
   F->InteriorCutoff = DBLSTRF (&InputStr);
   F->ExteriorCutoff = DBLSTRF (&InputStr);
   F->f0 = (FunctionCall_t *) GetTableValue;
   F->f1 = (FunctionCall_t *) GetTableDeriv;

   /*  Sanity Check  */

   if (F->NumCoeff <  2)
		{
		printf
         ("ERROR:  ");
      printf
         ("There must be at least 2 data points in table potenial.\n");
      CleanAfterError();
      }

   if (F->ExteriorCutoff <= F->InteriorCutoff)
      {
      printf ( "ERROR:  Exterior Cutoff (%e)",
         F->ExteriorCutoff);
		printf (" is less than or equal to interior cutoff (%e).\n",
			F->InteriorCutoff);
      CleanAfterError();
      }

   F->TableFactor =
      F->NumCoeff / ( F->ExteriorCutoff - F->InteriorCutoff );

   /*  Make permanant space for coefficients  */
   ALLOCATE (F->Coeff, double, 6*F->NumCoeff)

   /*  With data array in hand, call table interpolation routine  */
	InterpolateTable(F, a);

	}
#endif


void InterpolateTable (Function_t *F, double *a)
	{
	double *p;
	int i,j,k,l,m,n;

	/*  Condition temporary array according to boundary conditions  */
	switch (GET_INT_EXTRAP(F->BoundaryType))
		{
		case FUNC_EXTRAP_CONST_SLOPE:
			a[ -1] = 2.0*a[  0] - a[1];
			a[ -2] = 2.0*a[ -1] - a[0];
			break;
		case FUNC_EXTRAP_CONST_VALUE:
			a[ -1] = a[  0];
			a[ -2] = a[  0];
			break;
		case FUNC_EXTRAP_ZERO:
			a[ -1] = 0.0;
			a[ -2] = 0.0;
			break;
		default:
			/*  Assertion will print error if this is reached  */
			ASSERT(0)
		}

	n = F->NumCoeff+1;
   switch (GET_EXT_EXTRAP(F->BoundaryType))
      {
      case FUNC_EXTRAP_CONST_SLOPE:
         a[n  ] = 2.0*a[n-1] - a[n-2];
         a[n+1] = 2.0*a[n  ] - a[n-1];
         break;
      case FUNC_EXTRAP_CONST_VALUE:
         a[n  ] = a[n-1];
			a[n+1] = a[n-1];
         break;
      case FUNC_EXTRAP_ZERO:
         a[n  ] = 0.0;
         a[n+1] = 0.0;
         break;
      default:
			/*  Assertion will print error if this is reached  */
			ASSERT(0)
		}


	/*   INTERPOLATE DERIVATIVES  */
	p = F->Coeff;
   LOOP (k, F->NumCoeff)
      {
      i = k-2;
      j = k-1;
      l = k+1;
      m = k+2;
      n = k+3;
      *p++ = (                      a[k]                            )   ;
      *p++ = ( 2*a[i] -16*a[j]           + 16*a[l] - 2*a[m]         )/24;
      *p++ = (  -a[i] +16*a[j] - 30*a[k] + 16*a[l] -   a[m]         )/24;
      *p++ = (-9*a[i] +39*a[j] - 70*a[k] + 66*a[l] -33*a[m] + 7*a[n])/24;
      *p++ = (13*a[i] -64*a[j] +126*a[k] -124*a[l] +61*a[m] -12*a[n])/24;
      *p++ = (-5*a[i] +25*a[j] - 50*a[k] + 50*a[l] -25*a[m] + 5*a[n])/24;
      }

	/*  Calculate slope at endpoints  */
	i = 6*0;
	p = F->Coeff;
	F->InteriorValue = a[0];
	F->InteriorSlope = F->TableFactor* p[i+1];

	i = 6*F->NumCoeff;
	F->ExteriorValue = a[F->NumCoeff];
	F->ExteriorSlope = F->TableFactor*
		(p[i+1] + 2*p[i+2] + 3*p[i+3] + 4*p[i+4] + 5*p[i+5]);

   }


/*
Read Polynomial Function from input
*/
static void ReadPolyFunction
(
Function_t *F,
char       *InputStr,
LIST       *InList
)
   {
   char ReadBuffer[NBUFF];
   char *TempStr;
   BOOLEAN EndOfFile;
   int iread;
   int icoeff;

   F->Power1         = INTSTRF (&InputStr);
   F->Power2         = INTSTRF (&InputStr);
   F->x0             = DBLSTRF (&InputStr);
   F->ExteriorCutoff = DBLSTRF (&InputStr);
   F->InteriorCutoff = 0.0;
   F->NumCoeff       = F->Power2 - F->Power1 + 1;
   F->f0             = (FunctionCall_t *) GetPolyValue;
   F->f1             = (FunctionCall_t *) GetPolyDeriv;

   /*  Allocate space for coefficients  */
	if (F->Coeff!=NULL)
		FREE(F->Coeff)
   ALLOCATE (F->Coeff, double, 2*F->NumCoeff)

   /*  Read coefficients  */
   iread = 0;
   ReadBuffer[0] = '\0';
   TempStr          = ReadBuffer;
   EndOfFile      = FALSE;
   while (iread<F->NumCoeff && !EndOfFile)
      {

      /*  Read new line if current line is blank  */
      if (IsBlank(TempStr))
         {
         EndOfFile =
            (NULL==FGETS(ReadBuffer, NBUFF, InList));
         TempStr = ReadBuffer;
         }

      /*  Premature end-of-file, end program  */
      if (EndOfFile)
         {
         printf ("ERROR: File ends before coefficients were found\n");
         CleanAfterError();
         }

      F->Coeff[iread] = DBLSTRF (&TempStr);
      iread++;
      }

   /*  Calcualte coefficients for first derivative  */
   LOOP (icoeff, F->NumCoeff)
      {
      F->Coeff [icoeff + F->NumCoeff] =
         ( icoeff + F->Power1 ) * F->Coeff[icoeff];
      }

	/*  Calculate slope at endpoints  */
	F->InteriorValue = GetPolyValue(F, F->InteriorCutoff);
	F->InteriorSlope = GetPolyDeriv(F, F->InteriorCutoff);
	F->ExteriorValue = GetPolyValue(F, F->ExteriorCutoff);
	F->ExteriorSlope = GetPolyDeriv(F, F->ExteriorCutoff);

   }



static void ScaleTableFunction (Function_t *F, double Scale, int ScaleType)
   {
   int iseg;
   int ipow;
   int icoeff;

	/*   
	To scale function input, 
	only need to scale factor for converting argument units
   to table units (table nodes spaced a unit distance apart)
	*/
   if (FUNC_INPUT==ScaleType)
      {
      F->TableFactor /= Scale;
      }

   else if (FUNC_OUTPUT==ScaleType)
      {
      LOOP (icoeff, 6*F->NumCoeff)
         F->Coeff[icoeff] *= Scale;
      }
   }


static void ScalePolyFunction (Function_t *F, double Scale, int ScaleType)
   {
   int icoeff;

   /*
   Scale function input
   */
   if (FUNC_INPUT==ScaleType)
      {

      /*  Scale polynomial origin  */
      F->x0             *= Scale;

      /*  Scale DERIV=0 Coefficients  */
      LOOP (icoeff, F->NumCoeff)
         {
         F->Coeff[icoeff] /= pow(Scale, F->Power1+icoeff);
         }

      /*  Scale DERIV=1 Coefficients  */
      LOOP (icoeff, F->NumCoeff)
         {
         F->Coeff[icoeff+F->NumCoeff]
            /= pow(Scale, F->Power1+icoeff);
         }

      ASSERT (F->Power1+F->NumCoeff==F->Power2)
      }

   /*
   Scale function output
   */
   else if (FUNC_OUTPUT==ScaleType)
      {
      LOOP (icoeff, 2*F->NumCoeff)
         {
         F->Coeff[icoeff] *=Scale;
         }
      }

   /*
   ERROR
   */
   else
      {
      printf ("INTERNAL ERROR:  Invalid scale flag.\n");
      CleanAfterError();
      }
   }


static double GetTableValue (Function_t *F, double x)
   {
   double *p;
   double q;
   int    m;

   if (NULL==F->Coeff)
      return 0.0;

   if (x<F->InteriorCutoff || x>=F->ExteriorCutoff)
      return OutOfRangeValue (F, x, DERIV_ZERO);

   m  = (q = F->TableFactor * (x - F->InteriorCutoff) );
   q -= m;

   /*
   This should happen because out of range values are handled above
   */
   ASSERT(m<F->NumCoeff)

   /*  Point to m'th item (6 numbers per item)  */
   p = &(F->Coeff[6*m]);

   /*  Interpolate  */
   return p[0] + (p[1] + (p[2] + (p[3] + (p[4] + p[5]*q)*q)*q)*q)*q;

   }


static double GetTableDeriv (Function_t *F, double x)
   {
   int m;
   double q;
   double *p;

	if (F->Coeff==NULL)
      return 0.0;

   if (x<F->InteriorCutoff || x>=F->ExteriorCutoff)
      return OutOfRangeValue (F, x, DERIV_ONE);

   m  = (q = F->TableFactor * (x - F->InteriorCutoff) );
   q -= m;

   /*
   This should happen because out of range values are handled above
   */
   ASSERT(m<F->NumCoeff)

   /*  Point to m'th item (6 numbers per item)  */
   p = &(F->Coeff[6*m]);

   /*  Interpolate  */
   return
      F->TableFactor *
      ( p[1] + (2*p[2] + (3*p[3] + (4*p[4] + 5*p[5]*q)*q)*q)*q );
   }



static double GetPolyValue (Function_t *F, double x)
   {
   double dx;
   double *p;

   if (x<=F->InteriorCutoff || x>=F->ExteriorCutoff)
      return OutOfRangeValue (F, x, DERIV_ZERO);

   dx = x - F->x0;
   p  = F->Coeff;

   if (0==F->Power1)
      return poly (dx, F->NumCoeff-1, p);
   else if (1==F->Power1)
      return dx * poly (dx, F->NumCoeff-1, p);
   else if (2==F->Power1)
      return dx * dx * poly (dx, F->NumCoeff-1, p);
   else
      return pow (dx, F->Power1) * poly (dx, F->NumCoeff-1, p);
   }



static double GetPolyDeriv (Function_t *F, double x)
   {
   double dx;
   double *p;

   if (x<F->InteriorCutoff || x>F->ExteriorCutoff)
      return OutOfRangeValue (F, x, DERIV_ONE);

   dx = x - F->x0;
   p  = F->Coeff+F->NumCoeff;

   if (0==F->Power1)
      return poly (dx, F->NumCoeff-2, p+1);
   else if (1==F->Power1)
      return poly (dx, F->NumCoeff-1, p);
   else if (2==F->Power1)
      return  dx * poly (dx, F->NumCoeff-1, p);
   else
      return pow (dx, F->Power1-1) * poly (dx, F->NumCoeff-1, p);
   }



static double GetNullValue (Function_t *F, double x)
	{
	return 0.0;
	}



static double OutOfRangeValue (Function_t *F, double x, int Deriv)
   {
   double dx;


   /*
   Argument less than interior cutoff
   Do this first because POLY only has this kind of cutoff
   */
   if (x <= F->InteriorCutoff)
      {
      switch (GET_INT_BOUND(F->BoundaryType))
         {

         case FUNC_BOUND_CONST_SLOPE:
            if (DERIV_ONE==Deriv)
               {
               return F->InteriorSlope;
               }
            else if (DERIV_ZERO==Deriv)
               {
               dx  = x - F->InteriorCutoff;
               return F->InteriorValue + F->InteriorSlope * dx;
               }
				break;

			case FUNC_BOUND_CONST_VALUE:
				if (DERIV_ONE==Deriv)
					{
					return 0.0;
					}
				else
					{
					return F->InteriorValue;
					}

			case FUNC_BOUND_ZERO:
				return 0.0;

			case FUNC_BOUND_OUTOFRANGE:
				ERROR_PREFIX
				printf ("Function within interior cutoff.\n");
				CleanAfterError();
				break;

			}

      ASSERT(0)
      }

   /*  Argument exceeds exterior cutoff  */
   if (x >= F->ExteriorCutoff)
      {
      switch (GET_EXT_BOUND(F->BoundaryType))
         {
         case FUNC_BOUND_CONST_SLOPE:
            if (DERIV_ONE==Deriv)
               {
               return F->ExteriorSlope;
               }
            else if (DERIV_ZERO==Deriv)
               {
               dx  = x - F->ExteriorCutoff;
               return F->ExteriorValue + F->ExteriorSlope * dx;
               }
            break;
         case FUNC_BOUND_CONST_VALUE:
            if (DERIV_ONE==Deriv)
               {
               return F->ExteriorValue;
               }
            else
               {
               return 0.0;
               }
			case FUNC_BOUND_ZERO:
            return 0.0;
			case FUNC_BOUND_OUTOFRANGE:
				ERROR_PREFIX
				printf ("Function within exterior cutoff.\n");
            CleanAfterError();
            break;
         }
		ASSERT(0)
      }



   /*  Execution should never reach here  */
	ASSERT (0)

   return 0.0;
   }



double poly (double x, int n, double *p)
   {
   double ReturnValue = p[n];

   while (n>0)
      {
      n--;
      ReturnValue = ReturnValue * x + p[n];
      }

   return ReturnValue;
   }



/*
************************************************************************
Local Function for Stand-Alone Test Version
************************************************************************
*/
#ifdef LOCAL



BOOLEAN IsWhiteSpace (char InChar)
   {
   return
      InChar==' '  ||
      InChar=='\t' ||
      InChar=='\r' ||
      InChar=='\n';
   }



BOOLEAN IsBlank (char *String)
   {
   while (IsWhiteSpace(*String))
      String++;
   return (*String=='\0');
   }



char *strhed (char **String)
   {
   char *TokenStr;

   /*  Find first non-whitespace  */
   while (IsWhiteSpace(**String))
      {
      (*String)++;
      }
   TokenStr = *String;

   /*  Find Next Whitespace Character  */
   while (!IsWhiteSpace(**String) && **String!='\0')
      {
      (*String)++;
      }

   /*  If end of line, return  */
   if (**String=='\0')
      {
      return TokenStr;
      }

   /*  Mark as end of string, find next non-whitespace  */
   (**String)='\0';
   (*String)++;
   while (IsWhiteSpace(**String))
      (*String)++;

   return TokenStr;
   }


double dblstr (char **String)
   {
   char *TokenStr;

   TokenStr = strhed (String);
   return atof (TokenStr);
   }

double intstr (char **String)
   {
   return atoi (strhed(String));
   }

void CleanAfterError(void)
   {
   exit(0);
   }

BOOLEAN IsFirstToken(char *Token, char *String)
	{
	int Len;
	while (IsWhiteSpace(*String))
		String++;
	return 
		!STRNCMPI(Token,String,Len) && 
		(IsWhiteSpace(String[Len]) || String[Len]=='\0');
	}




/*
************************************************************************
Test Driver
************************************************************************
*/

int main ()
   {
   Function_t *F = NULL;
   int BoundaryType;
   char Buffer[1024];
   char *InStr  = NULL;
   char *TokStr = NULL;
   char *String = "POLY 2 2 4 4";
   double Cutoff;
   int i;
   double x;


   strcpy (Buffer, String);
   inlist = stdin;

	/*  Intialize Function  */
	F = FUNC_CreateNullFunction();
   F->BoundaryType =
      MAKE_BOUND_FLAG
         (
         FUNC_EXTRAP_CONST_SLOPE, FUNC_BOUND_OUTOFRANGE,
         FUNC_EXTRAP_ZERO,        FUNC_BOUND_ZERO
         );

   fgets (Buffer, 1024, stdin);
   InStr = Buffer;
   TokStr = strhed (&InStr);
   FUNC_ParseRead (F, InStr);

   Cutoff = FUNC_GetCutoff(F);

   LOOP (i, 101)
      {
      x = F->InteriorCutoff + 0.01*i*(F->ExteriorCutoff-F->InteriorCutoff);
      printf ("%le %le %le\n", x, F->f0(F,x), F->f1(F,x));
      }

/*  Write out values  */
#if 0
   printf ("x f(x): %le %le\n", 0.0, FUNC_GetValue(F,0.0));
   printf ("x f(x): %le %le\n", 0.5, FUNC_GetValue(F,0.5));
   printf ("x f(x): %le %le\n", 2.0, FUNC_GetValue(F,2.0));
   printf ("x f(x): %le %le\n", 4.0, FUNC_GetValue(F,4.0));

   printf ("x f0(x): %le %le\n", 0.0, F->f0(F,0.0));
   printf ("x f1(x): %le %le\n", 0.0, F->f1(F,0.0));
#endif

   return 0;
   }

#endif