csmbase.c

/**************************************************************/
/*                           IDCP                             */
/*            The Insertion Device Control Program            */
/*                 Author: Götz Pfeiffer                     */
/**************************************************************/

/*____________________________________________________________*/
/* project information:                                       */
/*____________________________________________________________*/
/* project: IDCP                                              */
/* module: CSM-BASE                                           */
/* version-number of module: 1.1                              */
/* author: Götz Pfeiffer                                     */
/* last modification date: 2011-08-16                         */
/* status: tested                                             */
/*____________________________________________________________*/

/*____________________________________________________________*/
/*                          comments                          */
/*____________________________________________________________*/

    /*@ETI------------------------------------------------*/
    /*     comments for the automatically generated       */
    /*  header-file, this does not concern csm.c but      */
    /*                      csmbase.h !                   */
    /*----------------------------------------------------*/

/*@EM("/@ This file was automatically generated, all manual changes\n")
  @EM("   to this file will be lost, edit csm.c instead of csm.h !@/\n\n") */

    /*@ETI------------------------------------------------*/
    /*			    History                       */
    /*----------------------------------------------------*/

/*
Version 0.9:
  Date: 6/16/97
  This is the first implementation csm.

Version 0.91:
  Date: 8/22/97
  A mechanism to accelerate the search of x-values in the table has been
  implemented. Now the interval-borders from the previous search are
  tested (and if it makes sense taken as borders for the new search).
  This should speed up the process if table_x_lookup_ind is called with
  x-values that do not differ too much from each other (in terms of the
  table-values). This version has not yet been tested.

Version 0.92:
  Date: 10/16/97
  The program (csmbase) has been splitted from the original csm and now
  contains only generic routines for implementing linear and table-based
  conversions. All IDCP specific conversion funtions are now placed
  in csm.c. Note that for this file (csmbase.c) version 0.92 is the first
  version!

Version 0.93:
  Date: 5/5/99
  csmbase was extended in order to support 2-dimensional functions (z=f(x,y))
  These functions are specified in form of a table :
	 y1   y2   y3   .   .   .
     x1  z11  z12  z13  .   .   .
     x2  z21  z22  z23  .   .   .
     .    .    .    .
     .    .    .    .
     .    .    .    .

Version 0.94:
  Date: 6/16/99
  The functions csm_read_table and csm_read_xytable were changed in order to
  work even if some lines of the data-file do not contain data (this is
  especially the case with empty lines).

Version 0.95:
  Date: 5/26/00
  An error in coordinate_lookup_index() was removed. This caused the function
  to return wrong parameters when applied with x==2.

Version 0.96:
  Date: 10/13/00
  Small changes. If a function is not yet defined (type CSM_NOTHING), the
  function returns 0.

  Date: 2002-03-05
  Errors in the header-file generation were fixed, some
  left-overs from previous debugging were removed

  Date: 2004-04-05
  in csm_read_table a superfluous parameter was removed

  Date: 2004-07-23
  a new function, csm_clear was added. This function releases all
  memory an csm_function occupies and sets its value to CSM_NOTHING

  Date: 2004-08-31
  a new function, csm_free was added. This function releases all
  memory and finally releases the memory of the csm_function object
  itself.

  Date: 2006-03-01
  bugfix: semDelete was added

  Date: 2006-03-01
  * usage of psem (portable semaphores) can now be selected
  * usage of DBG, errlogPrinf or printf can now be selected
    define USE_DBG and USE_ERRLOGPRINTF accordingly
  * prepared for Linux (undefine B_VXWORKS for this,
    set USE_DBG 0, USE_ERRLOGPRINTF 0 and USE_PSEM 1
  * memory leak error was fixed in csm_read_2d_table

  Date: 2006-03-02
  cosmetic changes

  Date: 2006-03-02
  changes in the embedded documentation

  Date: 2007-11-14
  * bugfix with respect to the on-hold function: The "last" field
    of the csm_Function structure has to exist for each type of
    function that can be applied to the structure, one for
    csm_x(), one for csm_y(), one for csm_z(), for csm_dx() and
    csm_dy(). With only one field, it may happen for example, 
    that when "on_hold" is used, csm_x() returns the last returned 
    value of csm_y(). 
  * bugfix with respect to multithreading use:
    the semaphore used in functions like csm_x(), csm_y() and so
    on was given just before a return-statement which did a read-access
    on the function structure. This lead to a race condition where
    it was possible that one call of csm_x() returned the value of 
    another call of csm_x() which was called by another task. Together
    with "last" field bug mentioned above, csm_x() could even return
    the value of the last csm_y() call, which is fatal...
  * str_empty_or_comment() now accepts lines where "#" is not 
    in the first column.  
  * the inline documentation was updated  

  Date: 2007-11-14
  * the usage of semaphores can no be disabled by a macro.
    This was needed for better portability.

  Date: 2011-08-16
  * The functions lookup_1d_functiontable and lookup_2d_functiontable now
    return NAN (not a number) in case the table doesn't contain any values. 
  * lookup_1d_functiontable now handles the case of a table with a single
    point (xp,yp) correctly. It then returns the associated yp for any value of
    x. Before this change, it returned NAN for x!=xp and yp for x==xp.
*/

    /*----------------------------------------------------*/
    /*		        General Comments                  */
    /*----------------------------------------------------*/

/*! \file csmbase.c
    \brief implement function table
    \author Götz Pfeiffer (goetz.pfeiffer\@helmholtz-berlin.de)
    \version 1.0

   This module implements one and two dimensional function tables.
   The tables can be specified by simple ascii files. The points
   defined in the file don't have to be in the same distance.
*/


/*@ITI________________________________________________________*/
/*                      general Defines                       */
/*@ET_________________________________________________________*/

/*@EM("\n#ifndef __CSMBASE_H\n#define __CSMBASE_H\n") */

#ifndef DOXYGEN_SKIP_THIS

#ifndef B_LINUX
/*! \internal \brief select compilation under vxWorks */
#define B_VXWORKS
#endif

/*! \internal \brief needed for POSIX compability */
#define _INCLUDE_POSIX_SOURCE

/*! \internal \brief use DBG module ? */
#define USE_DBG 0

/*! \internal \brief use epics errlogPrintf ? */
#define USE_ERRLOGPRINTF 1

/*! \internal \brief use semaphores at all ? */
#define USE_SEM 1

/*! \internal \brief use psem semaphore module ? */
#define USE_PSEM 0

/* the following macros affect the debug (dbg) module: */

#if USE_DBG

/*! \internal \brief compile DBG assertions */
#define DBG_ASRT_COMP 0

/*! \internal \brief compile DBG debug messages */
#define DBG_MSG_COMP 0

/*! \internal \brief compile DBG trace messages */
#define DBG_TRC_COMP  0

/*! \internal \brief use local switch-variables for the dbg module */
#define DBG_LOCAL_COMP 1

/*! \internal \brief use async. I/O with message passing in DBG module*/
#define DBG_ASYNC_COMP 1

/* the following defines are for sci-debugging, they do not influence
   any header-files but are placed here since they are usually changed
   together with the above macros.*/

/*! \internal \brief for debug-outputs, "stdout" means console */
#define CSMBASE_OUT_FILE "stdout"

/*! \internal \brief csm trace level

  \li level 6: dense traces for debugging
  \li level 5: enter/exit tracing
*/
#define CSMBASE_TRC_LEVEL 0

/*! \internal \brief flushmode, 0, 1 and 2 are allowed */
#define CSMBASE_FLUSHMODE 1

#else

#if !USE_ERRLOGPRINTF
#define errlogPrintf printf
#endif

/*! \internal \brief compability macro, needed when DBG is not used */
#define DBG_MSG_PRINTF2(f,x) errlogPrintf(f,x)

/*! \internal \brief compability macro, needed when DBG is not used */
#define DBG_MSG_PRINTF3(f,x,y) errlogPrintf(f,x,y)

/*! \internal \brief compability macro, needed when DBG is not used */
#define DBG_MSG_PRINTF4(f,x,y,z) errlogPrintf(f,x,y,z)

/*! \internal \brief compability macro, needed when DBG is not used */
#define DBG_MSG_PRINTF5(f,x,y,z,q) errlogPrintf(f,x,y,z,q)

#endif

#if !USE_SEM
#define SEM_TYPE int
#define SEMTAKE(x) 
#define SEMGIVE(x) 
#define SEMDELETE(x) 
#define SEMCREATE(x) 0 

#else

#if USE_PSEM
#define SEM_TYPE psem_mutex
#define SEMTAKE(x) psem_mutex_take(&(x))
#define SEMGIVE(x) psem_mutex_give(&(x))
#define SEMDELETE(x) psem_mutex_remove(&(x))
#define SEMCREATE(x) (!psem_mutex_create(&(x)))
#else
/*! \internal \brief compability macro for semaphores */
#define SEM_TYPE epicsMutexId
/*! \internal \brief compability macro for semaphores */
#define SEMTAKE(x) epicsMutexLock(x)
/*! \internal \brief compability macro for semaphores */
#define SEMGIVE(x) epicsMutexUnlock(x)
/*! \internal \brief compability macro for semaphores */
#define SEMDELETE(x) epicsMutexDestroy(x)
/*! \internal \brief compability macro for semaphores */
#define SEMCREATE(x) (NULL==(x= epicsMutexCreate()))
#endif

#endif

#endif

/*! \brief for type csm_bool */
#define CSM_TRUE 1

/*! \brief for type csm_bool */
#define CSM_FALSE 0

/*____________________________________________________________*/
/*			 Include-Files			      */
/*____________________________________________________________*/

#if USE_SEM

#if USE_PSEM
#include <psem.h>
#else
#include <epicsMutex.h>
#endif

#endif

#if USE_DBG
#include <dbg.h>
#endif

#if USE_ERRLOGPRINTF
#include <errlog.h> /* epics error printf */
#endif

#include <ctype.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

/* math.h for definition of NAN */
#include <math.h>

/*@ITI________________________________________________________*/
/*                         Defines                            */
/*@ET_________________________________________________________*/

      /*................................................*/
      /*@IL           the NAN number type               */
      /*................................................*/

#ifndef NAN
/*! \internal \brief define NAN if it is not yet defined */
#define NAN (-(0.0/0.0))
#endif

/*@ITI________________________________________________________*/
/*                          Types                             */
/*@ET_________________________________________________________*/

      /*................................................*/
      /*@IL        the csm_bool type (public)           */
      /*................................................*/

/*! \brief the boolean data type used in this module */

typedef int csm_bool; /*@IL*/

      /*................................................*/
      /*          the coordinate type (private)         */
      /*................................................*/

/*! \internal \brief a single coordinate value */
typedef struct
  {
    double value; /*!< the floating point value of the coordinate */
    int    index; /*!< the index-number within this list of coordinates */
  } csm_coordinate;

/*! \internal \brief a list of coordinates */
typedef struct
  { csm_coordinate *coordinate;  /*!< pointer to coordinate array */
    int            no_of_elements;
                                 /*!< number of elements in the coordinate
				      array */
    int            a_last;       /*!< last left index */
    int            b_last;       /*!< last right index */
    int		   ret_last;     /*!< last return-code of lookup function */
    double         x_last;       /*!< last x that was looked up */
  } csm_coordinates;

      /*................................................*/
      /*        the linear function type (private)      */
      /*................................................*/

/*! \internal \brief a linear function y = a + b*x */
typedef struct
  { double a;         /*!< a in y= a+b*x */
    double b;         /*!< b in y= a+b*x */
  } csm_linear_function;

      /*................................................*/
      /*           the 1d table type (private)          */
      /*................................................*/

/*! \internal \brief one-dimenstional function table

  This is a list of x-y pairs that define a one-dimensional
  function. The function can be inverted, if it is monotone.
*/
typedef struct
  { csm_coordinates x;    /*!< list of x-coordinates */
    csm_coordinates y;    /*!< list of y-coordinates */
  } csm_1d_functiontable;

      /*................................................*/
      /*          the 2d table type (private)           */
      /*................................................*/

/*! \internal \brief two-dimenstional function table

  This is a list of x-y pairs and corresponding z-values that
  define a two-dimensional function table. The function cannot be inverted.
*/
typedef struct
  { csm_coordinates x;   /*!< list of x-coordinates */
    csm_coordinates y;   /*!< list of y-coordinates */
    double          *z;  /*!< list of z-values (function values) */
  } csm_2d_functiontable;

      /*................................................*/
      /*           the function-type (private)          */
      /*................................................*/

#ifndef DOXYGEN_SKIP_THIS
/*! \internal \brief typedef-enum: types of functions known in this module */
typedef enum { CSM_NOTHING,  /*!< kind of NULL value */
               CSM_LINEAR,   /*!< linear y=a+b*x */
               CSM_1D_TABLE, /*!< one-dimensional function table */
	       CSM_2D_TABLE  /*!< two-dimensional function table */
	     } csm_func_type;
#endif

/*! \internal \brief compound type for functions */
struct csm_Function
  { SEM_TYPE semaphore;   /*!< semaphore to lock access to csm_function */
    double   last_x;      /*!< last x value that was calculated */
    double   last_y;      /*!< last y value that was calculated */
    double   last_dx;     /*!< last dx value that was calculated */
    double   last_dy;     /*!< last dy value that was calculated */
    double   last_z;      /*!< last z value that was calculated */
    csm_bool on_hold;     /*!< when TRUE, just return last */
    csm_func_type type;   /*!< the type of the function */

    union { csm_linear_function    lf;   /*!< linear function */
            csm_1d_functiontable   tf_1; /*!< 1d function table */
            csm_2d_functiontable   tf_2; /*!< 2d function table */
	  } f;            /*!< union that holds the actual data */

  };

      /*................................................*/
      /*@IL       the function-object (public)          */
      /*................................................*/

/*! \brief the abstract csm function object */
/*@IT*/
typedef struct csm_Function csm_function;
/*@ET*/

/*____________________________________________________________*/
/*			  Variables			      */
/*____________________________________________________________*/

      /*................................................*/
      /*           initialization (private)             */
      /*................................................*/

/*! \internal \brief initialization flag [static]

  This variable is used to remember, wether the module
  was already initialized.
*/
static csm_bool initialized= CSM_FALSE;

/*____________________________________________________________*/
/*                        Functions                           */
/*@ET_________________________________________________________*/


    /*----------------------------------------------------*/
    /*             debug - managment (private)            */
    /*----------------------------------------------------*/

#if USE_DBG
/*! \internal \brief initialize the DBG module [static]

  This internal function is called by csmbase_init
*/

static void csm_dbg_init(void)
  { static csm_bool csmbase_init= CSM_FALSE;

    if (csmbase_init)
      return;
    csmbase_init= CSM_TRUE;
    /* initialize the debug (dbg) module: use stdout for error-output */
    DBG_SET_OUT(CSMBASE_OUT_FILE);
    DBG_TRC_LEVEL= CSMBASE_TRC_LEVEL;
                           /* level 6: dense traces for debugging
                              level 5: enter/exit tracing
                              level 2: output less severe errors and warnings
                              level 1: output of severe errors
                              level 0: no output */
    DBG_FLUSHMODE(CSMBASE_FLUSHMODE);
  }

#endif

    /*----------------------------------------------------*/
    /*        utilities for csm-coordinates (private)     */
    /*----------------------------------------------------*/

        /*              initialization                */

/*! \internal \brief initialize a csm_coordinates structure [static]

  This function initializes a csm_coordinates structure. This function
  should be called immediately after a variable of the type
  csm_coordinates was created.
  \param c pointer to the csm_coordinates array
*/
static void init_coordinates(csm_coordinates *c)
  {
    c->a_last=-1;
    c->b_last=-1;
    c->x_last=-1;
    c->ret_last=-1;
    c->coordinate= NULL;
    c->no_of_elements=0;
  }

/*! \internal \brief initialize a csm_coordinates structure [static]

  This function allocates memory for a
  csm_coordinates structure.
  \param c pointer to the csm_coordinates structure
  \param elements number of elements in the structure
  \return returns CSM_FALSE in case of an error, CSM_TRUE otherwise
*/
static csm_bool alloc_coordinates(csm_coordinates *c, int elements)
  { int i;
    csm_coordinate *co;

    if (c->no_of_elements==0) /* 1st time memory allocation */
      c->coordinate= malloc(sizeof(csm_coordinate)*elements);
    else
      c->coordinate= realloc(c->coordinate,
                             sizeof(csm_coordinate)*elements);

    if (c->coordinate==NULL) /* allocation error */
      { DBG_MSG_PRINTF2("error in csm:alloc_coordinates line %d,\n"
                	"allocation failed!\n", __LINE__);
	init_coordinates(c);
	return(CSM_FALSE);
      };

    for(i=c->no_of_elements, co= (c->coordinate + c->no_of_elements);
        i<elements;
	i++, co++)
      { co->value=0;
	co->index=i;
      };

    c->no_of_elements= elements;
    return(CSM_TRUE);
  }

/*! \internal \brief resize a csm_coordinates structure [static]

  This function resizes a csm_coordinates structure. This means that
  the size of allocated memory is adapted to a new, different number
  of elements. Note that the existing data is not changed. If additional
  memory is allocated, it is initialized with zeroes.
  \param c pointer to the csm_coordinates structure
  \param elements new number of elements in the structure
  \return returns CSM_FALSE in case of an error, CSM_TRUE otherwise
*/
static csm_bool resize_coordinates(csm_coordinates *c, int elements)
  { if (elements==c->no_of_elements)
      return(CSM_TRUE);

    c->no_of_elements= elements;
    if (NULL== (c->coordinate= realloc(c->coordinate,
                                       sizeof(csm_coordinate)*elements)))
      { DBG_MSG_PRINTF2("error in csm:resize_coordinates line %d,\n" \
                        "realloc failed!\n", __LINE__);
        return(CSM_FALSE);
      };
    return(CSM_TRUE);
  }

/*! \internal
    \brief re-initialize a csm_coordinates structure [static]

  This re-initializes a csm_coordinates structure. This function
  is similar to init_coordinates, but already allocated memory is
  freed before the structure is filled with it's default values.
  \param c pointer to the csm_coordinates array
*/
static void reinit_coordinates(csm_coordinates *c)
  { if (c->no_of_elements==0)
      return;

    free(c->coordinate);
    init_coordinates(c);
  }

/*! \internal \brief initialize a matrix of double values [static]

  This function allocates an array of doubles in a way that it
  has enough room to hold all values of a rows*columns matrix.
  \param z address of pointer to array of double-numbers (output)
  \param rows number of rows
  \param columns number of columns
  \return returns CSM_FALSE in case of an error, CSM_TRUE otherwise
*/
static csm_bool init_matrix(double **z, int rows, int columns)
  { int i;
    double *zp;

    if (NULL== (zp= malloc((sizeof(double))*rows*columns)))
      { DBG_MSG_PRINTF2("error in csm:init_matrix line %d,\n" \
                        "malloc failed!\n", __LINE__);
        return(CSM_FALSE);
      };
    *z = zp;
    for (i= rows*columns; i>0; i--, *(zp++)=0);
    return(CSM_TRUE);
  }

        /*            x-compare for qsort             */

/*! \internal
    \brief compare two values within a csm_coordinate structure [static]

  This function compares two values within a cms_coordinate structure.
  It is needed in order to apply quicksort.
  \param t1 this is the pointer to the first value. It should be of the
            type csm_coordinate*.
  \param t2 this is the pointer to the second value. It should be of the
            type csm_coordinate*.
  \return returns -1 if *t1 < *t2, 0 if *t1 = *t2, 1 if *t1 > *t2
*/
static int coordinate_cmp(const void *t1, const void *t2)
  { double x1= ((csm_coordinate *)(t1))->value;
    double x2= ((csm_coordinate *)(t2))->value;

    if (x1 < x2)
      return(-1);
    if (x1 > x2)
      return( 1);
    return(0);
  }

        /*              coordinate-sort               */

/*! \internal \brief sort coordinates in a csm_coordinates structure [static]

  This function sorts the coordinates within a
  csm_coordinates structure according to their value field.
  \param coords pointer to the csm_coordinates structure
*/
static void coordinate_sort( csm_coordinates *coords)
  {
    qsort( coords->coordinate, coords->no_of_elements,
           sizeof(csm_coordinate), coordinate_cmp);
  }

/*! \internal \brief update backlinks from on csm_coordinates structure to another [static]

  Two or more csm_coordinates structures can be connected by the
  index fields of their elements. The index field of each element
  in the first structure is then the index of the corresponding
  element in the second structure and vice versa. If the order
  of the element in one of the structures was changed, the index
  field in all elements of the second structure have to be updated,
  which is what this function does. It is usually called after
  sorting the first structure with coordinate_sort.
  \param coords1 pointer to the first csm_coordinates structure which
         was re-ordered.
  \param coords2 pointer to the second csm_coordinates structure which
         has to be updated.
*/
static void coordinate_update_backlinks( csm_coordinates *coords1,
				         csm_coordinates *coords2)
  { int i,l;
    csm_coordinate *c,*d;

    l= coords1->no_of_elements;
    for (i=0, c= coords1->coordinate, d= coords2->coordinate; i<l; i++, c++)
      { d[c->index].index = i; };
  }

        /*      lookup an index in coordinates        */

/*! \internal \brief lookup nearest boundaries within a csm_coordinates structure [static]

  This searches for the nearest boundaries to a given x within a given
  csm_coordinates structure. It looks up the indices of two
  elements. For the first one,it is guaranteed that its value is
  smaller or equal to x and that no other element with a larger value
  has this property. For the second one, it is guaranteed that its value
  is bigger or equal to x and that no other element with a smaller value
  has this property. so
     coords->coordinate[a].value <= x <= coords->coordinate[b].value
  holds. If a or b cannot be found since x is at the border of the
  list of coordinates, they are set to -1.
  \param coords pointer to the first csm_coordinates structure
  \param x the value to find
  \param a pointer to the returned index of the smaller element
  \param b pointer to the returned index of the bigger element
  \return 2 if x was found, a is the corresponding index
          1 x was not found, but a and b define the closest interval around x
          0 x was not found at it is outside the closest interval a and b
            a or b may be -1
         -1 error
*/
static int coordinate_lookup_index(csm_coordinates *coords,double x,
				   int *a, int *b)
  /* 2: x was found, it is returned in a
     1: x was not found, but a and b define the closest interval around x
     0: x was not found at it is outside the closest interval a and b
    -1: error
  */
  /* the table should be x-sorted */
  { int i, ia, ib;
    double xt;
    csm_coordinate *c= coords->coordinate;

    if (coords->a_last==-1)
      { coords->a_last=0;
        coords->b_last=coords->no_of_elements-1;
      }
    else
      { if (x==coords->x_last)
          { *a= coords->a_last; *b= coords->b_last; return(coords->ret_last); };
      };

    *a= 0; *b=0;
    if (coords->no_of_elements<2)
      { if (coords->no_of_elements<1)
          {
            DBG_MSG_PRINTF2("error in csm:coordinate_lookup_index %d,\n" \
                            "table has less than 1 element!\n", __LINE__);
            return(-1);
           };
        coords->a_last= 0; coords->b_last= 0; coords->x_last= x;
	if (x== c[0].value)
	  return( coords->ret_last= 2 );
	return( coords->ret_last= 0 );
      };
    ia= 0;
    ib= coords->no_of_elements-1;

    if (c[ia].value>=x) /* left of area */
      { if (c[ia].value==x) /* exact match */
          { coords->a_last= *a= ia;
	    coords->b_last= *b= ia;
	    coords->x_last= x;
	    return( coords->ret_last=2 );
	  };
        coords->a_last= *a= ia;
        coords->b_last= *b= (ia+1);
	coords->x_last= x;
	return( coords->ret_last=0 );
      };
    if (c[ib].value<=x) /* right of area */
      { if (c[ib].value==x) /* exact match */
          { coords->a_last= *a= ib;
	    coords->b_last= *b= ib;
	    coords->x_last= x;
	    return( coords->ret_last=2 );
	  };
        coords->a_last= *a= (ib-1);
        coords->b_last= *b= ib;
        coords->x_last= x;
        return( coords->ret_last=0 );
      };

    /* search acceleration for consecutive values of x: */
    i= (coords->a_last)-1;
    i= (i<ia) ? ia : i;
    xt= c[i].value;
    if (xt==x)
      { coords->a_last= *a= i;
        coords->b_last= *b= i;
        coords->x_last= x;
	return( coords->ret_last=2 );
      };
    if (xt<x)
      ia= i;
    else
      ib= i;
    i= (coords->b_last)+1;
    i= (i>ib) ? ib : i;
    xt= c[i].value;
    if (xt==x)
      { coords->a_last= *a= i;
        coords->b_last= *b= i;
        coords->x_last= x;
	return( coords->ret_last=2 );
      };
    if (xt<x)
      ia= i;
    else
      ib= i;

    for (; (ib>ia+1); )
      { i= (ib+ia)/2;
        xt= c[i].value;
        if (xt<x)
          ia= i;
        else
          if (xt>x)
            ib= i;
          else  /* exact match */
            { coords->a_last= *a= i;
              coords->b_last= *b= i;
              coords->x_last= x;
              return( coords->ret_last=2 );
            };
      };
    coords->a_last= *a= ia;
    coords->b_last= *b= ib;
    coords->x_last= x;
    return( coords->ret_last=1 );
  }

    /*----------------------------------------------------*/
    /*     utilities for csm_linear_function (private)    */
    /*----------------------------------------------------*/

        /*               calculation                  */

/*! \internal \brief compute y from a given x [static]

  This function computes y when x is given for a linear function.
  A linear function is a
  \param p pointer to the linear function structure
  \param x the value of x
*/
static double linear_get_y(csm_linear_function *p, double x)
  { return( (p->a) + (p->b)*x ); }

/*! \internal \brief compute x from a given y [static]

  This function computes x when y is given for a linear function
  \param p pointer to the linear function structure
  \param y the value of y
*/
static double linear_get_x(csm_linear_function *p, double y)
  { return( (y - (p->a))/(p->b) ); }

/*! \internal \brief compute delta-y from a given delta-x [static]

  This function computes delta-y when delta-x is given for a linear function
  \param p pointer to the linear function structure
  \param x the value of x
*/
static double linear_delta_get_y(csm_linear_function *p, double x)
  { return( (p->b)*x ); }

/*! \internal \brief compute delta-x from a given delta-y [static]

  This function computes delta-x when delta-y is given for a linear function
  \param p pointer to the linear function structure
  \param y the value of y
*/
static double linear_delta_get_x(csm_linear_function *p, double y)
  { return( y/(p->b) ); }

    /*----------------------------------------------------*/
    /*    utilities for csm_1d_functiontable (private)    */
    /*----------------------------------------------------*/

        /*              initialization                */

/*! \internal \brief initialize a csm_1d_functiontable structure [static]

  This function initializes all members of the csm_1d_functiontable
  structure. This function should be called immediately after a
  variable of the type csm_1d_functiontable was created.
  \param ft pointer to the csm_1d_functiontable structure
*/
static void init_1d_functiontable(csm_1d_functiontable *ft)
  { init_coordinates(&(ft->x));
    init_coordinates(&(ft->y));
  }

/*! \internal \brief reinitialize a csm_1d_functiontable structure [static]

  This function re-initializes the csm_coordinates members of the
  csm_1d_functiontable structure. This is similar to init_1d_functiontable,
  but already allocated memory is freed before the structures is filled
  with their default values.
  \param ft pointer to the csm_1d_functiontable structure
*/
static void reinit_1d_functiontable(csm_1d_functiontable *ft)
  { reinit_coordinates(&(ft->x));
    reinit_coordinates(&(ft->y));
  }

        /*        lookup a value in the table         */

/*! \internal \brief lookup a value in a csm_1d_functiontable structure [static]

  This function looks up the function value y for a given x in a
  csm_1d_functiontable structure which is a one-dimensional function
  table. The two points x1 and x2 that are
  nearest to x are searched. A linear interpolation between the
  corresponding y1 and y2 values is done and this value is returned.
  If inverted lookup is wanted. the function looks up y1 and y2 for a
  given y and does a linear interpolation between x1 and x2.
  \param ft pointer to the csm_1d_functiontable structure
  \param x the value to look up
  \param invert if CSM_TRUE, an inverted lookup is done
  \return the value that fits best to x is returned.
*/
static double lookup_1d_functiontable(csm_1d_functiontable *ft, double x,
				      csm_bool invert)
  { int res,a,b;
    double xa,xb,ya,yb;
    csm_coordinate c;
    csm_coordinates *xcoords;
    csm_coordinates *ycoords;

    if (!invert)
      { xcoords= &(ft->x);
        ycoords= &(ft->y);
      }
    else
      { xcoords= &(ft->y);
        ycoords= &(ft->x);
      };

    res= coordinate_lookup_index(xcoords, x, &a, &b);


    if (res==-1) /* error */
      return(NAN);
    if ((res== 2)||(a==b)) /* exact match or table with just a single value */
      { 
        /* if the table has just a single value (only one [x,y] pair) we return
           the corresponding y value. */
        c= (xcoords->coordinate)[a];
        return( (ycoords->coordinate)[c.index].value );
      };

    c= (xcoords->coordinate)[a];
    xa= c.value;
    ya= (ycoords->coordinate)[c.index].value;

    c= (xcoords->coordinate)[b];
    xb= c.value;
    yb= (ycoords->coordinate)[c.index].value;

    /* y= ya + (x-xa)/(xb-xa)*(yb-ya) */
    return( ya + (x-xa)/(xb-xa) * (yb-ya) );
  }

    /*----------------------------------------------------*/
    /*    utilities for csm_2d_functiontable (private)    */
    /*----------------------------------------------------*/

        /*              initialization                */

/*! \internal \brief initialize a csm_2d_functiontable structure [static]

  This function initializes all members of the csm_2d_functiontable
  structure. This function should be called immediately after a
  variable of the type csm_2d_functiontable was created.
  \param ft pointer to the csm_2d_functiontable structure
*/
static void init_2d_functiontable(csm_2d_functiontable *ft)
  { init_coordinates(&(ft->x));
    init_coordinates(&(ft->y));
  }

/*! \internal \brief reinitialize a csm_2d_functiontable structure [static]

  This function re-initializes the csm_coordinates members of the
  csm_2d_functiontable structure. This is similar to init_2d_functiontable,
  but already allocated memory is freed before the structures is filled
  with their default values.
  \param ft pointer to the csm_2d_functiontable structure
*/
static void reinit_2d_functiontable(csm_2d_functiontable *ft)
  { reinit_coordinates(&(ft->x));
    reinit_coordinates(&(ft->y));
    if (ft->z!=NULL)
      free(ft->z);
    ft->z= NULL;
  }

        /*       lookup a value in the xy-table       */

/*! \internal \brief lookup a value in a csm_2d_functiontable structure [static]

  This function looks up the function value z for a given x and y in a
  csm_2d_functiontable structure which is a two-dimensional function
  table. The four points (x1,y1), (x1,y2), (x2,y1) and (x2,y2) that are
  closest to (x,y) are searched. A two-dimensional linear interpolation
  between four corresponding z values is done and this value is returned.
  \param ft pointer to the csm_2d_functiontable structure
  \param x x coordinate of the point to look up
  \param y y coordinate of the point to look up
  \return the value that fits best to (x,y) is returned.
*/
static double lookup_2d_functiontable(csm_2d_functiontable *ft,
				      double x, double y)
  { int res,xia,xib,yia,yib;
    double xa,xb,ya,yb;
    double zaa,zab,zba,zbb;
    double *z= ft->z;
    int no_of_columns;

    csm_coordinate cxa,cxb,cya,cyb;