cdfunc.c

一个很好的分子动力学程序

			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