5.4.2 How to build your own potential descriptors

Although this subject really is one that should be deferred to Chapter , we will now present a simple prototype ``definition'' for a potential descriptor in C and Fortran^5.12

    void inipotential (int *npar, double *par, char *name);
    void potential (int *ndim, double *pos, double *acc, double *pot, double *time);

    SUBROUTINE INIPOTENTIAL(NPAR, PAR, NAME)
    SUBROUTINE POTENTIAL(NDIM, POS, ACC, POT, TIME)

As you can and will see more of, a potential descriptor is in origin really a C or FORTRAN source code file, that needs two (FORTRAN) subroutines or (C) functions with the callable names inipotential and potential. Their arguments must conform to the specification given above. Because we do want to allow Fortran source code as well, all arguments are called by reference in C.

Programs which need a potential descriptor will automatically compile your source code (if needed) and load the object code into the program for usage. The repository of standard NEMO potential descriptors (as object files) lives in $NEMOOBJ/potential, and is automatically searched when the environment variable POTPATH is appropriately set. Note that the two subroutines themselves are not called directly by the user, but by a workhorse routine from the standard NEMO library. This hides much of the interface for the programmer. More details on this technique can be found in Chapter (still to come).

Below is a fully commented listing of the harmonic potential, as an example of such a potential descriptor given in the C language. It, and other potentials, can be found in source code form in the directory $NEMO/src/orbit/potential/data

/*
 * harmonic.c: procedures for initializing and calculating 
 *             the forces and potential of a harmonic potential
 */
#include <stdinc.h>                 /* formal NEMO include file */
local double omega = 0.0;       /* defined but not used in here */
local double h[3] = {1.0,1.0,1.0};        /* default parameters */
/*----------------------------------------------------------------
 * INIPOTENTIAL: initializes the potential.
 *      input: npar, the number of parameters
 *             par[] an array of npar parameters
 *      If npar=0 defaults are taken (remember to initialize them 
 *      as static (local) variables in this file)
 *----------------------------------------------------------------
 */
void inipotential (int *npar, double *par, string name)
{
    int i;

    if (*npar>0)
        omega = par[0];
    for (i=1; i<(*npar); i++)
       h[i-1] = sqr(par[i]);
    if (*npar > 4) 
        warning("Only 4 parameters used in Harmonic Potential");

    dprintf (1,"INI_POTENTIAL Harmonic Potential\n");
    dprintf (1,"  Parameters : Pattern Speed = %f\n",omega);
    dprintf (1,"  wx^2,wy^2,wz^2= %f %f %f\n",h[0],h[1],h[2]);
}
/*----------------------------------------------------------------
 *  POTENTIAL: the worker routine. Determines at any given point 
 *      (x,y,z) the forces and potential. 
 *      Note that this routine is good for 1, 2 as well as 3D
 *----------------------------------------------------------------
 */
void potential (int ndim,double *pos,double *acc,double *pot,double *time)
{
    int    i;
        
    *pot = 0.0;
    for (i=0; i<*ndim; i++) {
        (*pot) += h[i]*sqr(pos[i]);
        acc[i] = -h[i]*pos[i];
    }
    *pot *= 0.5;
}