home *** CD-ROM | disk | FTP | other *** search
/ Dr. CD ROM (Annual Premium Edition) / premium.zip / premium / IBMOS2_1 / LSQRFT13.ZIP / fitutil2.c < prev    next >
Text File  |  1993-07-10  |  5KB  |  219 lines
  1. #include <math.h>
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <string.h>
  5. #include "fit.h"
  6.  
  7. /* This file contains a bunch of function for allocating */
  8. /* vectors and matrices.  */
  9.  
  10. extern int debug;
  11. void myerror(char *);
  12.  
  13. double **dmatrix(nr,nc)
  14. int nr,nc;
  15. {
  16.     int i;
  17.     double **m;
  18.     char s[80];
  19.  
  20.     m=(double **) malloc((unsigned) (nr)*sizeof(double*));
  21.     if (!m) myerror("allocation failure 1 in dmatrix()");
  22.  
  23.     for(i=0;i<nr;i++) {
  24.         m[i]=(double *) malloc((unsigned) (nc)*sizeof(double));
  25.         sprintf(s,"allocation failure 2 in dmatrix() on column %d",i);
  26.         if (!m[i]) myerror(s);
  27.     }
  28.     return m;
  29. }
  30.  
  31. int **imatrix(nr,nc)
  32. int nr,nc;
  33. {
  34.     int i,**m;
  35.  
  36.     m=(int **)malloc((unsigned) (nr)*sizeof(int*));
  37.     if (!m) myerror("allocation failure 1 in imatrix()");
  38.  
  39.     for(i=0;i<nr;i++) {
  40.         m[i]=(int *)malloc((unsigned) (nc)*sizeof(int));
  41.         if (!m[i]) myerror("allocation failure 2 in imatrix()");
  42.     }
  43.     return m;
  44. }
  45.  
  46. void free_dmatrix(m,nr,nc)
  47. double **m;
  48. int nr,nc;
  49. {
  50.     int i;
  51.  
  52.     for(i=nr-1;i>=0;i--) free((char*) (m[i]));
  53.     free((char*) (m));
  54. }
  55.  
  56. void free_imatrix(m,nr,nc)
  57. int **m;
  58. int nr,nc;
  59. {
  60.     int i;
  61.  
  62.     for(i=nr-1;i>=0;i--) free((char*) m[i]);
  63.     free((char*) (m));
  64. }
  65.  
  66. void myerror(s)
  67. char s[80];
  68. {
  69. printf("%s\n", s);
  70. }
  71.  
  72. double *dvector(int n){
  73. return (double *)malloc((unsigned) (n)*sizeof(double));
  74. }
  75.  
  76. int *ivector(int n){
  77. return (int *)malloc((unsigned) (n)*sizeof(int));
  78. }
  79.  
  80.  
  81. /* parse parses command line arguments */
  82.  
  83. int parse(command, cmdarg)
  84. char *command, cmdarg[NUM_ARGS][30];
  85. {
  86. int i=0, j=0, k, toomany=0;
  87.  
  88. while( command[i] != '\x0' ){
  89.     if(j >= NUM_ARGS){
  90.         toomany = 1;
  91.         break;
  92.     }
  93.     if( command[i] == ' ' ){
  94.         i++;
  95.         k = 0;
  96.         while(command[i] == ' ') i++;
  97.         while( command[i] != ' ' && k < 30 && command[i] != '\x0'){
  98.             if(command[i] != ' '){
  99.                 cmdarg[j][k] = command[i];
  100.                 k++;
  101.             }
  102.         i++;
  103.        }
  104.     cmdarg[j][k] = '\x0';
  105.     j++;
  106.     }
  107. if(j==0)i++;
  108. }
  109. if(toomany){ 
  110.     printf("Warning: there is a maximum of %d command arguments\n",NUM_ARGS);
  111.     printf("Change NUM_ARGS in fit.h and recompile\n");
  112. }
  113. if(k == 0) j--;
  114. if(debug){
  115.     printf("#args: %d\n",j);
  116.     for( i = 0; i < j; i++) printf("%s\n", cmdarg[i]);
  117. }
  118. return j;
  119. }
  120.  
  121. /* calc_yfit calculates the value of the */
  122. /* fitting function for all the values of */
  123. /* the independent variable */
  124.  
  125. int calc_yfit(func, data, order, num_indep, a, ndata, ma, chisqr)
  126. int (*func)();
  127. double **data;
  128. struct data_order order;
  129. int num_indep;
  130. double a[];
  131. int ndata, ma;
  132. double *chisqr;
  133. {
  134. int i, j;
  135. int failed;
  136. double *dyda;
  137. double *dydx;
  138. int *fita;
  139. double *y;
  140. double *x;
  141. int *dydx_flag;
  142. double sig2i;
  143.  
  144. y = data[order.yfit];
  145.  
  146. dyda = dvector(ma);
  147. fita = ivector(ma);
  148. x = dvector(num_indep);
  149. dydx = dvector(num_indep);
  150. dydx_flag = ivector(num_indep);
  151.  
  152. for(j = 0; j < num_indep; j++) dydx_flag[j] = 0;
  153. for(j = 0; j < ma; j++) fita[j] = 0;
  154.  
  155. *chisqr = 0;
  156. for(i=0; i < ndata; i++){
  157.     for(j = 0; j < num_indep; j++){
  158.         x[j] = data[order.x[j]][i];
  159.         if(debug > 1)
  160.             printf("i: %d order.x[%d]: %d x[%d]: %g\n", i, j, order.x[j], j, x[j]);
  161.     }
  162.     if(debug > 1) printf("about to call (*func)()\n");
  163.     failed = (*func)(x,a,&y[i],dyda,ma,0,fita,dydx_flag,dydx,data[order.y][i]);
  164.     if(failed) break;
  165.     sig2i = data[order.sig][i]*data[order.sig][i];
  166.     if(sig2i > 1e-30)
  167.         *chisqr += (y[i] - data[order.y][i])*(y[i] - data[order.y][i])/sig2i;
  168. }
  169. free(dyda);
  170. free(fita);
  171. free(x);
  172. free(dydx);
  173. free(dydx_flag);
  174. return failed;
  175. }
  176.  
  177.  
  178. /* help lists comands */
  179. int help(char topic[12], int maxlines){
  180. char inbuf[80];
  181. FILE *stream;
  182. int i;
  183. char buf[20];
  184. if(strcmp(topic,"") == 0)strcpy(topic,"COMMANDS");
  185.  
  186. stream = NULL;
  187. if(getenv("FITHELP") == NULL && (stream = fopen("fit.hlp","rt")) == NULL){
  188.     printf("help file not found, set environment variable FITHELP to point to fit.hlp\n");
  189.     return 1;
  190. }
  191.  
  192. if( (stream == NULL) && (stream = fopen(getenv("FITHELP"),"rt")) == NULL){
  193.     printf("help file not found, set environment variable FITHELP to point to fit.hlp\n");
  194.     return 1;
  195. }
  196.  
  197. while(fgets(inbuf,80,stream) != NULL){
  198.     if( strncmp(inbuf, topic, strlen(topic)) == 0){
  199.         i = 0;
  200.         printf("%s",inbuf);
  201.         while( strncmp( fgets(inbuf,80,stream), "END",3) ){
  202.             printf("%s",inbuf);
  203.             i++;
  204.             if(i == 20){ i = 0; gets(buf); }
  205.             if( i == maxlines) return 0;
  206.         }
  207.     }
  208. }
  209. fclose(stream);
  210. return 0;
  211. }
  212.  
  213. int nonzero(double *array, int ndata){
  214. int i;
  215. for(i = 0; i < ndata; i++) if(array[i] < 1e-60) return 0;
  216. return 1;
  217. }
  218.  
  219.