csmbase.c 57.8 KB
Newer Older
franksen's avatar
franksen committed
1
2
3
/**************************************************************/
/*                           IDCP                             */
/*            The Insertion Device Control Program            */
4
/*                 Author: Götz Pfeiffer                     */
franksen's avatar
franksen committed
5
6
7
8
9
10
11
/**************************************************************/

/*____________________________________________________________*/
/* project information:                                       */
/*____________________________________________________________*/
/* project: IDCP                                              */
/* module: CSM-BASE                                           */
pfeiffer's avatar
pfeiffer committed
12
/* version-number of module: 1.1                              */
13
/* author: Götz Pfeiffer                                     */
14
/* last modification date: 2011-08-16                         */
pfeiffer's avatar
pfeiffer committed
15
/* status: tested                                             */
franksen's avatar
franksen committed
16
17
18
19
20
21
22
23
24
25
26
27
/*____________________________________________________________*/

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

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

28
/*@EM("/@ This file was automatically generated, all manual changes\n")
franksen's avatar
franksen committed
29
30
31
32
33
  @EM("   to this file will be lost, edit csm.c instead of csm.h !@/\n\n") */

    /*@ETI------------------------------------------------*/
    /*			    History                       */
    /*----------------------------------------------------*/
34
35

/*
franksen's avatar
franksen committed
36
37
Version 0.9:
  Date: 6/16/97
38
  This is the first implementation csm.
franksen's avatar
franksen committed
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

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
70
71
  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
franksen's avatar
franksen committed
72
73
74
75
76
77
78
79
80
81
82
  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.
83

franksen's avatar
franksen committed
84
  Date: 2002-03-05
85
  Errors in the header-file generation were fixed, some
franksen's avatar
franksen committed
86
  left-overs from previous debugging were removed
87

88
89
  Date: 2004-04-05
  in csm_read_table a superfluous parameter was removed
90
91
92

  Date: 2004-07-23
  a new function, csm_clear was added. This function releases all
93
94
95
96
97
  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
98
  itself.
pfeiffer's avatar
pfeiffer committed
99
100
101

  Date: 2006-03-01
  bugfix: semDelete was added
102
103
104
105
106
107
108
109

  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
pfeiffer's avatar
pfeiffer committed
110
111
112

  Date: 2006-03-02
  cosmetic changes
113
114
115

  Date: 2006-03-02
  changes in the embedded documentation
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132

  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...
133
134
  * str_empty_or_comment() now accepts lines where "#" is not 
    in the first column.  
135
  * the inline documentation was updated  
136
137
138
139

  Date: 2007-11-14
  * the usage of semaphores can no be disabled by a macro.
    This was needed for better portability.
140
141
142
143
144
145
146

  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.
franksen's avatar
franksen committed
147
148
149
150
151
*/

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

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

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

franksen's avatar
franksen committed
163
164
165
166
167
168
169

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

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

170
171
#ifndef DOXYGEN_SKIP_THIS

172
#ifndef B_LINUX
173
/*! \internal \brief select compilation under vxWorks */
174
175
176
#define B_VXWORKS
#endif

177
/*! \internal \brief needed for POSIX compability */
franksen's avatar
franksen committed
178
179
#define _INCLUDE_POSIX_SOURCE

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

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

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

189
/*! \internal \brief use psem semaphore module ? */
190
#define USE_PSEM 0
191
192
193
194
195

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

#if USE_DBG

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

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

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

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

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

franksen's avatar
franksen committed
211
/* the following defines are for sci-debugging, they do not influence
212
   any header-files but are placed here since they are usually changed
franksen's avatar
franksen committed
213
   together with the above macros.*/
214

215
/*! \internal \brief for debug-outputs, "stdout" means console */
franksen's avatar
franksen committed
216
#define CSMBASE_OUT_FILE "stdout"
217

218
/*! \internal \brief csm trace level
219
220
221
222

  \li level 6: dense traces for debugging
  \li level 5: enter/exit tracing
*/
franksen's avatar
franksen committed
223
#define CSMBASE_TRC_LEVEL 0
224

225
/*! \internal \brief flushmode, 0, 1 and 2 are allowed */
franksen's avatar
franksen committed
226
227
#define CSMBASE_FLUSHMODE 1

228
229
#else

230
231
232
233
#if !USE_ERRLOGPRINTF
#define errlogPrintf printf
#endif

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

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

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

pfeiffer's avatar
pfeiffer committed
243
244
245
/*! \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)

246
247
#endif

248
249
250
251
252
253
254
255
256
#if !USE_SEM
#define SEM_TYPE int
#define SEMTAKE(x) 
#define SEMGIVE(x) 
#define SEMDELETE(x) 
#define SEMCREATE(x) 0 

#else

257
258
259
260
261
262
263
264
#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 */
265
#define SEM_TYPE epicsMutexId
266
/*! \internal \brief compability macro for semaphores */
267
#define SEMTAKE(x) epicsMutexLock(x)
268
/*! \internal \brief compability macro for semaphores */
269
#define SEMGIVE(x) epicsMutexUnlock(x)
270
/*! \internal \brief compability macro for semaphores */
271
#define SEMDELETE(x) epicsMutexDestroy(x)
272
/*! \internal \brief compability macro for semaphores */
273
#define SEMCREATE(x) (NULL==(x= epicsMutexCreate()))
274
275
#endif

276
277
#endif

278
279
#endif

280
281
282
283
284
/*! \brief for type csm_bool */
#define CSM_TRUE 1

/*! \brief for type csm_bool */
#define CSM_FALSE 0
franksen's avatar
franksen committed
285
286
287
288
289

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

290
#if USE_SEM
291

292
293
294
#if USE_PSEM
#include <psem.h>
#else
295
#include <epicsMutex.h>
296
#endif
franksen's avatar
franksen committed
297

298
#endif
299

300
#if USE_DBG
301
#include <dbg.h>
302
303
304
#endif

#if USE_ERRLOGPRINTF
305
#include <errlog.h> /* epics error printf */
306
#endif
franksen's avatar
franksen committed
307

pfeiffer's avatar
pfeiffer committed
308
#include <ctype.h>
franksen's avatar
franksen committed
309
#include <stdlib.h>
310
#include <stdio.h>
franksen's avatar
franksen committed
311
312
#include <string.h>

313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
/* 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

franksen's avatar
franksen committed
329
330
331
332
333
/*@ITI________________________________________________________*/
/*                          Types                             */
/*@ET_________________________________________________________*/

      /*................................................*/
334
      /*@IL        the csm_bool type (public)           */
franksen's avatar
franksen committed
335
336
      /*................................................*/

337
/*! \brief the boolean data type used in this module */
franksen's avatar
franksen committed
338

339
typedef int csm_bool; /*@IL*/
340

341
342
343
      /*................................................*/
      /*          the coordinate type (private)         */
      /*................................................*/
344

345
/*! \internal \brief a single coordinate value */
franksen's avatar
franksen committed
346
typedef struct
347
  {
348
349
    double value; /*!< the floating point value of the coordinate */
    int    index; /*!< the index-number within this list of coordinates */
franksen's avatar
franksen committed
350
  } csm_coordinate;
351
352

/*! \internal \brief a list of coordinates */
franksen's avatar
franksen committed
353
typedef struct
354
  { csm_coordinate *coordinate;  /*!< pointer to coordinate array */
franksen's avatar
franksen committed
355
    int            no_of_elements;
356
                                 /*!< number of elements in the coordinate
357
358
359
360
361
				      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 */
362
  } csm_coordinates;
franksen's avatar
franksen committed
363

364
365
      /*................................................*/
      /*        the linear function type (private)      */
franksen's avatar
franksen committed
366
367
      /*................................................*/

368
/*! \internal \brief a linear function y = a + b*x */
franksen's avatar
franksen committed
369
typedef struct
370
371
  { double a;         /*!< a in y= a+b*x */
    double b;         /*!< b in y= a+b*x */
372
373
374
375
376
377
  } csm_linear_function;

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

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

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

388
389
390
391
      /*................................................*/
      /*          the 2d table type (private)           */
      /*................................................*/

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

394
  This is a list of x-y pairs and corresponding z-values that
395
  define a two-dimensional function table. The function cannot be inverted.
396
*/
franksen's avatar
franksen committed
397
typedef struct
398
399
400
401
402
  { 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;

403
404
      /*................................................*/
      /*           the function-type (private)          */
franksen's avatar
franksen committed
405
406
      /*................................................*/

407
#ifndef DOXYGEN_SKIP_THIS
408
409
410
411
412
/*! \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 */
413
	     } csm_func_type;
414
#endif
415

416
/*! \internal \brief compound type for functions */
417
struct csm_Function
418
  { SEM_TYPE semaphore;   /*!< semaphore to lock access to csm_function */
419
420
421
422
423
    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 */
424
    csm_bool on_hold;     /*!< when TRUE, just return last */
425
    csm_func_type type;   /*!< the type of the function */
426

427
428
429
430
    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 */
431
432

  };
franksen's avatar
franksen committed
433

434
435
      /*................................................*/
      /*@IL       the function-object (public)          */
franksen's avatar
franksen committed
436
437
      /*................................................*/

438
/*! \brief the abstract csm function object */
franksen's avatar
franksen committed
439
/*@IT*/
franksen's avatar
typo    
franksen committed
440
typedef struct csm_Function csm_function;
franksen's avatar
franksen committed
441
/*@ET*/
franksen's avatar
franksen committed
442

443
444
445
/*____________________________________________________________*/
/*			  Variables			      */
/*____________________________________________________________*/
446

447
448
449
450
      /*................................................*/
      /*           initialization (private)             */
      /*................................................*/

451
452
453
454
455
/*! \internal \brief initialization flag [static]

  This variable is used to remember, wether the module
  was already initialized.
*/
456
static csm_bool initialized= CSM_FALSE;
franksen's avatar
franksen committed
457
458
459
460
461
462
463

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


    /*----------------------------------------------------*/
464
    /*             debug - managment (private)            */
franksen's avatar
franksen committed
465
466
    /*----------------------------------------------------*/

467
468
#if USE_DBG
/*! \internal \brief initialize the DBG module [static]
franksen's avatar
franksen committed
469

470
471
  This internal function is called by csmbase_init
*/
franksen's avatar
franksen committed
472

473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
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);
  }
489

490
#endif
franksen's avatar
franksen committed
491
492

    /*----------------------------------------------------*/
493
    /*        utilities for csm-coordinates (private)     */
franksen's avatar
franksen committed
494
495
496
497
    /*----------------------------------------------------*/

        /*              initialization                */

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

500
501
  This function initializes a csm_coordinates structure. This function
  should be called immediately after a variable of the type
502
503
504
505
  csm_coordinates was created.
  \param c pointer to the csm_coordinates array
*/
static void init_coordinates(csm_coordinates *c)
506
  {
507
508
509
510
511
512
513
514
515
516
    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]

517
518
  This function allocates memory for a
  csm_coordinates structure.
519
520
521
522
523
  \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)
franksen's avatar
franksen committed
524
525
  { int i;
    csm_coordinate *co;
526

pfeiffer's avatar
pfeiffer committed
527
    if (c->no_of_elements==0) /* 1st time memory allocation */
528
      c->coordinate= malloc(sizeof(csm_coordinate)*elements);
529
    else
530
531
      c->coordinate= realloc(c->coordinate,
                             sizeof(csm_coordinate)*elements);
532

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

540
541
    for(i=c->no_of_elements, co= (c->coordinate + c->no_of_elements);
        i<elements;
542
543
544
	i++, co++)
      { co->value=0;
	co->index=i;
franksen's avatar
franksen committed
545
      };
546

547
548
    c->no_of_elements= elements;
    return(CSM_TRUE);
549
  }
franksen's avatar
franksen committed
550

551
552
553
554
555
/*! \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
556
  memory is allocated, it is initialized with zeroes.
557
558
559
560
561
  \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)
franksen's avatar
franksen committed
562
  { if (elements==c->no_of_elements)
563
      return(CSM_TRUE);
564

franksen's avatar
franksen committed
565
566
567
    c->no_of_elements= elements;
    if (NULL== (c->coordinate= realloc(c->coordinate,
                                       sizeof(csm_coordinate)*elements)))
568
      { DBG_MSG_PRINTF2("error in csm:resize_coordinates line %d,\n" \
569
570
                        "realloc failed!\n", __LINE__);
        return(CSM_FALSE);
franksen's avatar
franksen committed
571
      };
572
    return(CSM_TRUE);
573
  }
574

575
/*! \internal
576
577
578
579
    \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
580
  freed before the structure is filled with it's default values.
581
582
583
584
585
  \param c pointer to the csm_coordinates array
*/
static void reinit_coordinates(csm_coordinates *c)
  { if (c->no_of_elements==0)
      return;
586

587
588
    free(c->coordinate);
    init_coordinates(c);
589
  }
franksen's avatar
franksen committed
590

591
592
593
594
595
/*! \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)
596
597
  \param rows number of rows
  \param columns number of columns
598
599
600
  \return returns CSM_FALSE in case of an error, CSM_TRUE otherwise
*/
static csm_bool init_matrix(double **z, int rows, int columns)
franksen's avatar
franksen committed
601
602
  { int i;
    double *zp;
603

franksen's avatar
franksen committed
604
    if (NULL== (zp= malloc((sizeof(double))*rows*columns)))
605
      { DBG_MSG_PRINTF2("error in csm:init_matrix line %d,\n" \
606
607
                        "malloc failed!\n", __LINE__);
        return(CSM_FALSE);
franksen's avatar
franksen committed
608
609
610
      };
    *z = zp;
    for (i= rows*columns; i>0; i--, *(zp++)=0);
611
    return(CSM_TRUE);
612
  }
franksen's avatar
franksen committed
613
614
615

        /*            x-compare for qsort             */

616
/*! \internal
617
    \brief compare two values within a csm_coordinate structure [static]
franksen's avatar
franksen committed
618

619
620
621
622
623
624
625
626
  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
*/
franksen's avatar
franksen committed
627
628
629
static int coordinate_cmp(const void *t1, const void *t2)
  { double x1= ((csm_coordinate *)(t1))->value;
    double x2= ((csm_coordinate *)(t2))->value;
630

franksen's avatar
franksen committed
631
632
633
634
635
636
637
638
639
    if (x1 < x2)
      return(-1);
    if (x1 > x2)
      return( 1);
    return(0);
  }

        /*              coordinate-sort               */

640
641
642
643
644
645
/*! \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
*/
franksen's avatar
franksen committed
646
647
static void coordinate_sort( csm_coordinates *coords)
  {
648
    qsort( coords->coordinate, coords->no_of_elements,
franksen's avatar
franksen committed
649
650
651
           sizeof(csm_coordinate), coordinate_cmp);
  }

652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
/*! \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.
*/
franksen's avatar
franksen committed
667
static void coordinate_update_backlinks( csm_coordinates *coords1,
668
				         csm_coordinates *coords2)
franksen's avatar
franksen committed
669
670
  { int i,l;
    csm_coordinate *c,*d;
671
672

    l= coords1->no_of_elements;
franksen's avatar
franksen committed
673
    for (i=0, c= coords1->coordinate, d= coords2->coordinate; i<l; i++, c++)
674
      { d[c->index].index = i; };
franksen's avatar
franksen committed
675
676
677
678
  }

        /*      lookup an index in coordinates        */

679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
/*! \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
*/
701
static int coordinate_lookup_index(csm_coordinates *coords,double x,
702
				   int *a, int *b)
franksen's avatar
franksen committed
703
704
705
706
707
708
709
710
711
  /* 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;
712

713
714
    if (coords->a_last==-1)
      { coords->a_last=0;
franksen's avatar
franksen committed
715
        coords->b_last=coords->no_of_elements-1;
716
      }
franksen's avatar
franksen committed
717
718
719
720
    else
      { if (x==coords->x_last)
          { *a= coords->a_last; *b= coords->b_last; return(coords->ret_last); };
      };
721

franksen's avatar
franksen committed
722
723
724
    *a= 0; *b=0;
    if (coords->no_of_elements<2)
      { if (coords->no_of_elements<1)
725
          {
726
            DBG_MSG_PRINTF2("error in csm:coordinate_lookup_index %d,\n" \
727
                            "table has less than 1 element!\n", __LINE__);
franksen's avatar
franksen committed
728
729
            return(-1);
           };
730
        coords->a_last= 0; coords->b_last= 0; coords->x_last= x;
franksen's avatar
franksen committed
731
	if (x== c[0].value)
732
	  return( coords->ret_last= 2 );
franksen's avatar
franksen committed
733
734
735
736
	return( coords->ret_last= 0 );
      };
    ia= 0;
    ib= coords->no_of_elements-1;
737

franksen's avatar
franksen committed
738
739
    if (c[ia].value>=x) /* left of area */
      { if (c[ia].value==x) /* exact match */
740
741
          { coords->a_last= *a= ia;
	    coords->b_last= *b= ia;
franksen's avatar
franksen committed
742
743
	    coords->x_last= x;
	    return( coords->ret_last=2 );
744
745
746
	  };
        coords->a_last= *a= ia;
        coords->b_last= *b= (ia+1);
franksen's avatar
franksen committed
747
748
749
750
751
	coords->x_last= x;
	return( coords->ret_last=0 );
      };
    if (c[ib].value<=x) /* right of area */
      { if (c[ib].value==x) /* exact match */
752
753
          { coords->a_last= *a= ib;
	    coords->b_last= *b= ib;
franksen's avatar
franksen committed
754
	    coords->x_last= x;
755
756
757
758
	    return( coords->ret_last=2 );
	  };
        coords->a_last= *a= (ib-1);
        coords->b_last= *b= ib;
franksen's avatar
franksen committed
759
        coords->x_last= x;
760
        return( coords->ret_last=0 );
franksen's avatar
franksen committed
761
762
763
764
765
766
767
      };

    /* search acceleration for consecutive values of x: */
    i= (coords->a_last)-1;
    i= (i<ia) ? ia : i;
    xt= c[i].value;
    if (xt==x)
768
769
      { coords->a_last= *a= i;
        coords->b_last= *b= i;
franksen's avatar
franksen committed
770
        coords->x_last= x;
771
	return( coords->ret_last=2 );
franksen's avatar
franksen committed
772
773
774
775
776
777
778
779
780
      };
    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)
781
782
      { coords->a_last= *a= i;
        coords->b_last= *b= i;
franksen's avatar
franksen committed
783
        coords->x_last= x;
784
	return( coords->ret_last=2 );
franksen's avatar
franksen committed
785
786
787
788
789
      };
    if (xt<x)
      ia= i;
    else
      ib= i;
790

franksen's avatar
franksen committed
791
792
793
794
795
796
797
798
799
    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 */
800
801
            { coords->a_last= *a= i;
              coords->b_last= *b= i;
franksen's avatar
franksen committed
802
              coords->x_last= x;
803
              return( coords->ret_last=2 );
franksen's avatar
franksen committed
804
805
            };
      };
806
807
    coords->a_last= *a= ia;
    coords->b_last= *b= ib;
franksen's avatar
franksen committed
808
809
    coords->x_last= x;
    return( coords->ret_last=1 );
810
811
  }

812
813
814
815
816
    /*----------------------------------------------------*/
    /*     utilities for csm_linear_function (private)    */
    /*----------------------------------------------------*/

        /*               calculation                  */
franksen's avatar
franksen committed
817

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

  This function computes y when x is given for a linear function.
821
  A linear function is a
822
  \param p pointer to the linear function structure
823
  \param x the value of x
824
825
826
827
828
829
830
831
*/
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
832
  \param y the value of y
833
834
835
836
837
838
839
840
*/
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
841
  \param x the value of x
842
843
844
845
846
847
848
849
*/
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
850
  \param y the value of y
851
852
853
854
855
856
857
858
859
860
*/
static double linear_delta_get_x(csm_linear_function *p, double y)
  { return( y/(p->b) ); }

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

        /*              initialization                */

861
862
863
864
865
866
867
/*! \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
*/
868
869
870
871
872
static void init_1d_functiontable(csm_1d_functiontable *ft)
  { init_coordinates(&(ft->x));
    init_coordinates(&(ft->y));
  }

873
874
875
876
877
878
879
880
/*! \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
*/
881
882
883
884
static void reinit_1d_functiontable(csm_1d_functiontable *ft)
  { reinit_coordinates(&(ft->x));
    reinit_coordinates(&(ft->y));
  }
885

886
        /*        lookup a value in the table         */
franksen's avatar
franksen committed
887

888
889
890
891
892
893
894
895
896
897
898
899
900
901
/*! \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.
*/
franksen's avatar
franksen committed
902
static double lookup_1d_functiontable(csm_1d_functiontable *ft, double x,
903
				      csm_bool invert)
franksen's avatar
franksen committed
904
905
906
  { int res,a,b;
    double xa,xb,ya,yb;
    csm_coordinate c;
907
    csm_coordinates *xcoords;
franksen's avatar
franksen committed
908
909
910
    csm_coordinates *ycoords;

    if (!invert)
911
      { xcoords= &(ft->x);
franksen's avatar
franksen committed
912
913
914
        ycoords= &(ft->y);
      }
    else
915
      { xcoords= &(ft->y);
franksen's avatar
franksen committed
916
917
        ycoords= &(ft->x);
      };
918

franksen's avatar
franksen committed
919
920
921
922
    res= coordinate_lookup_index(xcoords, x, &a, &b);


    if (res==-1) /* error */
923
924
925
926
927
928
      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];
franksen's avatar
franksen committed
929
        return( (ycoords->coordinate)[c.index].value );
930
931
      };

franksen's avatar
franksen committed
932
933
934
    c= (xcoords->coordinate)[a];
    xa= c.value;
    ya= (ycoords->coordinate)[c.index].value;
935

franksen's avatar
franksen committed
936
937
938
    c= (xcoords->coordinate)[b];
    xb= c.value;
    yb= (ycoords->coordinate)[c.index].value;
939

franksen's avatar
franksen committed
940
941
942
943
    /* y= ya + (x-xa)/(xb-xa)*(yb-ya) */
    return( ya + (x-xa)/(xb-xa) * (yb-ya) );
  }

944
945
946
    /*----------------------------------------------------*/
    /*    utilities for csm_2d_functiontable (private)    */
    /*----------------------------------------------------*/
franksen's avatar
franksen committed
947

948
949
        /*              initialization                */

950
951
952
953
954
955
956
/*! \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
*/
957
958
959
960
961
static void init_2d_functiontable(csm_2d_functiontable *ft)
  { init_coordinates(&(ft->x));
    init_coordinates(&(ft->y));
  }

962
963
964
965
966
967
968
969
/*! \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
*/
970
971
972
973
974
975
976
977
978
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       */
franksen's avatar
franksen committed
979

980
981
982
983
984
985
986
987
988
989
990
991
/*! \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.
*/
992
static double lookup_2d_functiontable(csm_2d_functiontable *ft,
993
				      double x, double y)
franksen's avatar
franksen committed
994
995
996
997
998
999
1000
  { 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;
For faster browsing, not all history is shown. View entire blame