csmbase.c

/*! \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.
*/

/*____________________________________________________________*/
/*                      general Defines                       */
/*____________________________________________________________*/


#ifndef DOXYGEN_SKIP_THIS

/*! \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>

#include "csmbase.h"

/*____________________________________________________________*/
/*                         Defines                            */
/*____________________________________________________________*/

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

/*____________________________________________________________*/
/*                          Types                             */
/*____________________________________________________________*/

      /*................................................*/
      /*          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 */
  } 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 */
    double x_last;         /*!< last x value looked up */
    double y_last;         /*!< last y value looked up */
    csm_bool value_cached; /*!< 1 if x_last, y_last valid */
  } 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) */
    double x_last;         /*!< last x value looked up */
    double y_last;         /*!< last y value looked up */
    double z_last;         /*!< last z value looked up */
    csm_bool value_cached; /*!< 1 if x_last, y_last valid */
  } 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 */

  };

/*____________________________________________________________*/
/*			  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                           */
/*____________________________________________________________*/


    /*----------------------------------------------------*/
    /*             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->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
*/

#define ret(idx_a, idx_b, ret_code) \
  *a= idx_a;\
  *b= idx_b;\
  coords->a_last= idx_a;\
  coords->b_last= idx_b;\
  return(ret_code)
  
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;
    double x_boundary;
    int last_idx= coords->no_of_elements-1;
    int a_last, b_last;
    int tries;

    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__);
            ret(0,0,-1);
           };
        coords->a_last= 0; coords->b_last= 0; 
	if (x== c[0].value)
          {
            ret(0,0,2);
          }
        else
          {
            ret(0,0,0);
          }
      };

    a_last= coords->a_last;
    b_last= coords->b_last;
    if (a_last<0)
      a_last= 0;
    if (b_last<0)
      b_last= last_idx;

    for(tries=1; tries<=3; tries++)
      {
        switch(tries)
          {
            case 1:
              break;
            case 2:
              b_last= a_last;
              a_last= a_last-1;
              break;
            case 3:
              b_last= a_last;
              a_last= 0;
              break;
          }

        x_boundary= c[a_last].value;
        if (x>x_boundary)
          break;
        if (x==x_boundary)
          {
            ret(a_last,a_last,2);
          }
        if (a_last==0)
          { /* left of area */
            ret(0,1, 0);
          }
      }

    for(tries=1; tries<=3; tries++)
      {
        switch(tries)
          {
            case 1:
              break;
            case 2:
              a_last= b_last;
              b_last= b_last+1;
              break;
            case 3:
              a_last= b_last;
              b_last= last_idx;
              break;
          }

        x_boundary= c[b_last].value;
        if (x<x_boundary)
          break;
        if (x==x_boundary)
          {
            ret(b_last,b_last,2);
          }
        if (b_last==last_idx)
          { /* right of area */
            ret(last_idx-1,last_idx, 0);
          }
      }

    ia= a_last;
    ib= b_last;
    /* necessary condition here:
       c[ia].value<= x <= c[ib].value
     */
    for (; (ib>ia+1); )
      { i= (ib+ia)/2;
        xt= c[i].value;
        if (xt<x)
          {
            ia= i;
            continue;
          }
        if (xt>x)
          {
            ib= i;
            continue;
          }
        /* exact match */
        ret(i,i,2);
      };
    ret(ia,ib,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));
    ft->value_cached= CSM_FALSE;
  }

/*! \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));
    ft->value_cached= CSM_FALSE;
  }

        /*        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)
      { 
        if ((ft->value_cached) && (ft->x_last==x))
          return(ft->y_last);
        xcoords= &(ft->x);
        ycoords= &(ft->y);
      }
    else
      { 
        if ((ft->value_cached) && (ft->y_last==x))
          return(ft->x_last);
        xcoords= &(ft->y);
        ycoords= &(ft->x);
      };

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


    if (res==-1) /* error */
      {
        ft->value_cached= CSM_FALSE;
        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];
        ft->value_cached= CSM_TRUE;
        if (!invert)
          {
            ft->x_last= x;
            return(ft->y_last= (ycoords->coordinate)[c.index].value);
          }
        else
          {
            ft->y_last= x;
            return(ft->x_last= (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;

    ft->value_cached= CSM_TRUE;

    if (!invert)
      {
        ft->x_last= x;
        /* y= ya + (x-xa)/(xb-xa)*(yb-ya) */
        return(ft->y_last= ya + (x-xa)/(xb-xa) * (yb-ya));
      }
    else
      {
        ft->y_last= x;
        /* y= ya + (x-xa)/(xb-xa)*(yb-ya) */
        return(ft->x_last= 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));
    ft->value_cached= CSM_FALSE;
  }

/*! \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;
    ft->value_cached= CSM_FALSE;
  }

        /*       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;
    csm_coordinates *xcoords= &(ft->x);
    csm_coordinates *ycoords= &(ft->y);

    if (ft->value_cached)
      {
        if ((ft->x_last==x) && (ft->y_last==y))
          return(ft->z_last);
      }

    no_of_columns= ycoords->no_of_elements;

    res= coordinate_lookup_index(xcoords, x, &xia, &xib);
    if (res==-1) /* error */
      {
        ft->value_cached= CSM_FALSE;
        return(NAN);
      }

    res= coordinate_lookup_index(ycoords, y, &yia, &yib);
    if (res==-1) /* error */
      {
        ft->value_cached= CSM_FALSE;
        return(NAN);
      }

    /* res==2 with exact match is currently not treated separately */

    cxa= (xcoords->coordinate)[xia];
    cxb= (xcoords->coordinate)[xib];
    cya= (ycoords->coordinate)[yia];
    cyb= (ycoords->coordinate)[yib];

    xa= cxa.value; /* with exact match, xa==xb or ya==yb is possible */
    xb= cxb.value;
    ya= cya.value;
    yb= cyb.value;

    zaa= z[ cxa.index * no_of_columns + cya.index ];
    zab= z[ cxa.index * no_of_columns + cyb.index ];
    zba= z[ cxb.index * no_of_columns + cya.index ];
    zbb= z[ cxb.index * no_of_columns + cyb.index ];

#if 0
printf("x:%f y:%f xa:%f xb:%f ya:%f yb:%f\n",x,y,xa,xb,ya,yb);
printf("zaa:%f zab:%f zba:%f zbb:%f\n",zaa,zab,zba,zbb);
#endif
#if 1
    { double alpha = (xb==xa) ? 0 : (x-xa)/(xb-xa);
      double beta  = (yb==ya) ? 0 : (y-ya)/(yb-ya);
      if (alpha==0)               /* xb==xa ? */
        { 
          ft->value_cached= CSM_TRUE;
          ft->x_last= x;
          ft->y_last= y;
          if (beta==0)            /* yb==ya ? */
	    return(ft->z_last= zaa);
	  return(ft->z_last= zaa+ beta *(zab-zaa) );
	};
      if (beta==0)                /* yb==ya ? */
        {
          ft->value_cached= CSM_TRUE;
          ft->x_last= x;
          ft->y_last= y;
          return(ft->z_last= zaa+ alpha*(zba-zaa) );
        }

#if 0
printf(" ret: %f\n",
       zaa + beta*(zab-zaa)+alpha*(zba-zaa+beta*(zbb-zba-zab+zaa)));
#endif
      ft->value_cached= CSM_TRUE;
      ft->x_last= x;
      ft->y_last= y;
      return(ft->z_last= 
              zaa + beta*(zab-zaa)+alpha*(zba-zaa+beta*(zbb-zba-zab+zaa)) );
    };
#endif

#if 0
    { /* alternative according to formula given by J.Bahrdt */
      /* doesn't work when yb==ya or xb==xa */
      double delta= (yb-ya)*(xb-xa);
      double yb_  = yb - y;
      double xb_  = xb - x;
      double y_a  = y  - ya;
      double x_a  = x  - xa;
      return((zaa*yb_*xb_ + zab*y_a*xb_ + zba*yb_*x_a + zbb*y_a*x_a)/delta);
#endif

 }

    /*----------------------------------------------------*/
    /*              generic functions (public)            */
    /*----------------------------------------------------*/

double csm_x(csm_function *func, double y)
  { double ret;

    SEMTAKE(func->semaphore);

    if (func->on_hold)
      { ret= func->last_x;
        SEMGIVE(func->semaphore);
        return(ret);
      };

    switch (func->type)
      { case CSM_NOTHING:
               func->last_x= 0;
	       break;
	case CSM_LINEAR:
	       func->last_x= linear_get_x(&(func->f.lf), y);
	       break;
	case CSM_1D_TABLE:
	       func->last_x= lookup_1d_functiontable(&(func->f.tf_1),y,
	       					   CSM_TRUE);
	       break;
        default:
	       func->last_x=0;