home *** CD-ROM | disk | FTP | other *** search
/ Dr. CD ROM (Annual Premium Edition) / premium.zip / premium / IBMOS2_1 / LSQRFT13.ZIP / fitcmds2.c < prev    next >
Text File  |  1993-07-10  |  7KB  |  263 lines
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include "fit.h"
  5.  
  6. extern int debug;
  7.  
  8. /* function to get data and return number of data points */
  9. int get_data(data, command, num_indep, inbuf, order,filename,maxrows, maxcols)
  10. double **data;
  11. char *command;     /* command which was typed at the fit> prompt */
  12. int num_indep;
  13. char *inbuf;           /* inbuf is really output buffer */
  14. struct data_order *order;
  15. char *filename;
  16. int maxrows;
  17. int maxcols;
  18. {
  19. int i, j, jmax, ndata;
  20. char str[20][30];
  21. FILE *instream;
  22. int iopt1;
  23. char  inbuf1[512];
  24. iopt1 = 20;
  25.  
  26. sscanf(command,"%s %s %d", inbuf, filename, &iopt1);
  27.  
  28. if(( instream=fopen(filename,"rt") ) == NULL) return 0;
  29. i = 0;
  30. while( fgets(inbuf1,512,instream) != NULL && i < maxrows){
  31.     jmax=sscanf(inbuf1,"%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s",
  32.             str[0],str[1],str[2],str[3],str[4],str[5],str[6],str[7],str[8],str[9],
  33.             str[10],str[11],str[12],str[13],str[14],str[15],str[16],str[17],str[18],str[19]);
  34.     for(j=0; j < jmax && j <= iopt1 && j < maxcols-3; j++){
  35.         data[j][i] = atof(str[j]);
  36.     }
  37.     i++;
  38. }
  39. fclose(instream);
  40. ndata = i;
  41.  
  42. /* assign first columns to be independent variables */
  43. for(i=0; i < num_indep; i++)order->x[i] = i;
  44.  
  45. /* assign next column to be dependent variables */
  46. order->y = num_indep;
  47.  
  48. /* if there is a column num_indep+1, assign it to be error in y */
  49. if(jmax > num_indep) order->isig = num_indep + 1;
  50.  
  51. /* assign a data column to be all ones: no weighting */
  52. for(i = 0; i < ndata; i++)    data[jmax][i] = 1;
  53. order->nsig = jmax+0;
  54.  
  55. /* assign a data column to be used for statistical weighting */
  56. order->ssig= jmax+2;
  57.  
  58. /* assign a data column to be used for yfit */
  59. order->yfit= jmax+1;
  60.  
  61. if(jmax == 2 && num_indep == 1) order->sig = 2;
  62. if(jmax == 3 && num_indep == 2) order->sig = 3;
  63.  
  64. /* Assigning data columns for different weightings */
  65. /* takes up extra memory, but is faster and simpler */
  66. /* especially in the case of statistical weighting: */
  67. /* we don't have to take the sqrt of y at every point */
  68. /* for every iteration.  It is done once when statistical */
  69. /* weighting is selected */
  70.  
  71. return ndata;
  72. }
  73.  
  74. /* This function initializes the parameters */
  75. void ip(command, inbuf, func, a, ma)
  76. char *command, *inbuf;
  77. int (*func)();
  78. double *a;
  79. int ma;
  80. {
  81. int i;
  82. char cmdargs[NUM_ARGS][30];
  83.  
  84. if(func == NULL)
  85.     printf("ip failed, you must select a function first\n");
  86. else if(parse(command,cmdargs) == ma){
  87.     for( i = 0; i < ma; i++){
  88.         /* a * in place of a number leaves the parameter unchanged */
  89.         if(strcmp(cmdargs[i],"*") != 0) a[i] = atof(cmdargs[i]);
  90.     }
  91. }
  92. else
  93.     printf("ip failed, there should be %d parameter values entered\n",ma);
  94. }
  95.  
  96. /* This function changes one of the parameters */
  97. void cp(command, inbuf, func, a, ma)
  98. char *command, *inbuf;
  99. int (*func)();
  100. double *a;
  101. int ma;
  102. {
  103. int i;
  104. char cmdargs[NUM_ARGS][30];
  105.  
  106. if(func == NULL)
  107.     printf("cp failed, you must select a function first\n");
  108. else if(parse(command,cmdargs) == 2 && atoi(cmdargs[0]) > 0 
  109.         && atoi(cmdargs[0]) < ma)
  110.             a[atoi(cmdargs[0])] = atof(cmdargs[1]);
  111. else
  112.     help("cp",1);
  113. }
  114.  
  115. /* sp selects the parameters to vary */
  116. int sp(command, inbuf, func, lista, ma)
  117. char *command, *inbuf;
  118. int (*func)();
  119. int *lista, ma;
  120. {
  121. int i, mfit, max=0, min = 0;
  122. char cmdargs[NUM_ARGS][30];
  123.  
  124. if(func == NULL){
  125.     printf("sp failed, you must select a function first\n");
  126.     return 0;
  127. }
  128. else{
  129.     mfit = parse(command,cmdargs);
  130.     if(debug) printf("mfit: %d ma: %d\n", mfit, ma);
  131.     /* if mfit > ma, default to fitting all parameters */
  132.     if(ma < mfit){
  133.         mfit = ma;
  134.         for(i=0; i < ma; i++) lista[i] = i;
  135.         printf("sp failed, too many parameters selected\n");
  136.         if(debug)printf("mfit: %d ma: %d\n", mfit, ma);
  137.     }
  138.     else{
  139.         /* assign lista[i] to command arguments */
  140.         for( i = 0; i < mfit; i++) lista[i] = atoi(cmdargs[i]);
  141.         /* if any lista[i]'s are too big or too small, Default to all parameters */
  142.         for(i = 0; i < mfit; i++){
  143.             if( lista[i] > max) max = lista[i];
  144.             if( lista[i] < min) min = lista[i];
  145.         }
  146.         if(max > ma-1 || min < 0){
  147.             printf("sp failed, parameter out of range, all parameters fitted\n");
  148.             mfit = ma;
  149.             for(i=0; i < ma; i++) lista[i] = i;
  150.         }
  151.     }
  152.     }
  153. return mfit;
  154. }
  155.  
  156. int make_data(int (*func)(), int num_indep, double *a, int ma, char *command, char *inbuf){
  157.  
  158. char cmdargs[NUM_ARGS][30];
  159. char filename[80];
  160. FILE *stream;
  161. double *x, *xmin, *xmax, *xstep, y;
  162. int i,j, *fita;
  163. double  *dydx, *dyda;
  164. int *dydx_flag;
  165. int failed;
  166.  
  167. /* parse and check number of args entered on command line */
  168. if( (parse(command,cmdargs)) != 1 + 3*num_indep || func == NULL){
  169.     help("md",1);
  170.     return 1;
  171. }
  172.  
  173. if( (stream = fopen(cmdargs[0],"wt")) == NULL){
  174.     printf("md error, cannot open file %s\n", cmdargs[0]);
  175.     return 1;
  176. }
  177.  
  178. if( num_indep > 2 ){
  179.     printf("md error, md not implemented for more then 2 independent parameters\n");
  180.     return 1;
  181. }
  182.  
  183. /* need to send all these to the function */
  184. dydx = dvector(num_indep);
  185. dydx_flag = ivector(num_indep);
  186. fita = ivector(ma);
  187. dyda = dvector(ma);
  188.  
  189. /* we only need function value, not derivatives */
  190. for(j = 0; j < num_indep; j++) dydx_flag[j] = 0;
  191. for(j = 0; j < ma; j++) fita[j] = 0;
  192.  
  193. x = dvector(num_indep);
  194. xmin = dvector(num_indep);
  195. xmax = dvector(num_indep);
  196. xstep = dvector(num_indep);
  197.  
  198. /* assign command line arguments */
  199. for( j = 0; j < num_indep; j++){
  200.     x[j] = xmin[j] = atof(cmdargs[3*j+1]);
  201.     xmax[j] = atof(cmdargs[3*j+2]);
  202.     xstep[j] = atof(cmdargs[3*j+3]);
  203.     if(debug)printf("xmin: %g xmax: %g xstep: %g\n", xmin[j], xmax[j], xstep[j]);
  204. }
  205.  
  206. if(num_indep == 1){
  207.     while( x[0] <= xmax[0] ){
  208.         failed = (*func)(x,a,&y,dyda,ma,0,fita,dydx_flag,dydx);
  209.         if(failed){
  210.             free(dydx);
  211.             free(dydx_flag);
  212.             free(fita);
  213.             free(dyda);
  214.             free(x);
  215.             free(xmin);
  216.             free(xmax);
  217.             free(xstep);
  218.         return 1;
  219.         }
  220.         fprintf(stream,"%g %g\n", x[0], y);
  221.         x[0] += xstep[0];
  222.     }
  223. }
  224.  
  225. if(num_indep == 2){
  226.     while( x[0] <= xmax[0] ){
  227.         x[1] = xmin[1];
  228.         while( x[1] <= xmax[1] ){
  229.             failed = (*func)(x,a,&y,dyda,ma,0,fita,dydx_flag,dydx);
  230.             if(failed){
  231.                 free(dydx);
  232.                 free(dydx_flag);
  233.                 free(fita);
  234.                 free(dyda);
  235.                 free(x);
  236.                 free(xmin);
  237.                 free(xmax);
  238.                 free(xstep);
  239.                 return 1;
  240.             }
  241.             fprintf(stream,"%g %g %g \n", x[0], x[1], y);
  242.             if(debug)printf("%g %g %g\n", x[0], x[1], y);
  243.             x[1] += xstep[1];
  244.         }
  245.         x[0] += xstep[0];
  246.     }
  247. }
  248.  
  249. fclose(stream);
  250. strcpy(inbuf,command);
  251.  
  252. free(dydx);
  253. free(dydx_flag);
  254. free(fita);
  255. free(dyda);
  256. free(x);
  257. free(xmin);
  258. free(xmax);
  259. free(xstep);
  260. return 0;
  261.  
  262. }
  263.