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