util.cpp 3.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117
  1. /* util.c */
  2. /* Functions that are used for utility throughout the routines used for the C/S
  3. calculation. Most of these functions are from Numerical Recipes in C by Press. */
  4. #include "defs.h"
  5. #define NR_END 1
  6. #define FREE_ARG char*
  7. extern char errstr[200]; // error string that all routines have access to
  8. void nrerror(char error_text[])
  9. /* Numerical Recipes standard error handler */
  10. {
  11. fprintf(stderr,"Numerical Recipes run-time error...\n");
  12. fprintf(stderr,"%s\n",error_text);
  13. fprintf(stderr,"...now exiting to system...\n");
  14. exit(1);
  15. }
  16. float *fvector(int nl, int nh)
  17. /* allocate a float vector with subscript range v[nl..nh] */
  18. {
  19. float *v;
  20. v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
  21. if (!v) nrerror("allocation failure in fvector()");
  22. return v-nl+NR_END;
  23. }
  24. void free_fvector(float *v, int nl, int nh)
  25. /* free a float vector allocated with fvector() */
  26. {
  27. free((FREE_ARG) (v+nl-NR_END));
  28. }
  29. float **fmatrix(int nrl, int nrh, int ncl, int nch)
  30. /* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
  31. {
  32. int i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  33. float **m;
  34. /* allocate pointers to rows */
  35. m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
  36. if (!m) nrerror("allocation failure 1 in matrix()");
  37. m += NR_END;
  38. m -= nrl;
  39. /* allocate rows and set pointers to them */
  40. m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
  41. if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
  42. m[nrl] += NR_END;
  43. m[nrl] -= ncl;
  44. for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
  45. /* return pointer to array of pointers to rows */
  46. return m;
  47. }
  48. void free_fmatrix(float **m, int nrl, int nrh, int ncl, int nch)
  49. /* free a float matrix allocated by fmatrix() */
  50. {
  51. free((FREE_ARG) (m[nrl]+ncl-NR_END));
  52. free((FREE_ARG) (m+nrl-NR_END));
  53. }
  54. int copy_grid_geometry(FLOAT_GRID *grid_old, FLOAT_GRID *grid_new)
  55. // Copies geometric grid information from grid_old to grid_new
  56. {
  57. grid_new->start.x = grid_old->start.x;
  58. grid_new->start.y = grid_old->start.y;
  59. grid_new->start.z = grid_old->start.z;
  60. grid_new->inc.x = grid_old->inc.x;
  61. grid_new->inc.y = grid_old->inc.y;
  62. grid_new->inc.z = grid_old->inc.z;
  63. grid_new->x_count = grid_old->x_count;
  64. grid_new->y_count = grid_old->y_count;
  65. grid_new->z_count = grid_old->z_count;
  66. return(SUCCESS);
  67. }
  68. int binSearch(float *a, float searchnum, int M)
  69. /* Accepts a float array of data of length M ordered lowest to highest and a number called searchnum.
  70. Returns the index of the first element of the array, a, that is less than the searchnum.
  71. If the searchnum is less than a[0], then -1 is returned, and if the searchnum is greater
  72. than or equal to M, then M is returned. */
  73. {
  74. int found, mid, top, bottom;
  75. bottom = 0;
  76. top = M-1;
  77. found = 0; // flag that is set to 1 once the proper index is found
  78. // Ensure that the search parameter lies inside boundaries
  79. if(searchnum >= a[top])
  80. return(M);
  81. if(searchnum <= a[bottom])
  82. return(-1);
  83. while(!found)
  84. {
  85. mid = (top + bottom) / 2;
  86. if(searchnum == a[mid])
  87. found = 1;
  88. else
  89. if(searchnum < a[mid])
  90. top = mid - 1;
  91. else
  92. if(searchnum > a[mid + 1])
  93. bottom = mid + 1;
  94. else
  95. found = 1;
  96. }
  97. return(mid);
  98. }