home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 8 / CDASC08.ISO / NEWS / 2299 / GNUFIT / FIT2.C < prev    next >
Text File  |  1993-10-07  |  22KB  |  714 lines
  1. /*** fit2.c last modified 9 Jul 1993 ***/
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <string.h>
  5. #include <math.h>
  6. #include "fit.h"
  7.  
  8. /* GLOBAL DATA */
  9.  
  10. int grflag = 1;       /* do we graph data and fit? */
  11. int gnuopen = 0;      /* is there a stream to gnuplot open */
  12. FILE *pstream;        /* stream to gnuplot */
  13. double *xmin, *xmax;  /* min and max independent variables used in windowing */
  14. int wiflag = 0;       /* windowing flag */
  15. int veflag = 1;       /* verbosity flag */
  16. int debug = 0;        /* debugging flag */
  17. char GNUPLOT[60];     /* gnuplot command line */
  18.  
  19. main(int argc, char **argv){
  20.  
  21. double *a;            /* array of parameters for fit */
  22. double **covar;       /* covariant matrix for fit */
  23. double chisq;         /* squared error for fit */
  24. double **data;        /* matrix which holds data */
  25. int i, j, k;          /* indices for loops */
  26. int ndata=0;          /* number of data points */
  27. int ma;               /* number of parameters */
  28. int *lista;           /* array which has the list of parameters to vary */
  29. int mfit;             /* number of parameters being varied */
  30. int  itmax=20;        /* number of iterations */
  31. int iopt1, iopt2;     /* integers used in command options */
  32. char file[30];        /* name of file for reading or writing parameters*/
  33. char inbuf[80]="";    /* buffers used for input and output of strings */
  34. char buf[25];  
  35. char fname[20];        /* function name */
  36. char command[COMMAND_SIZE]="";   /* what was typed at command line */
  37. char cmdargs[NUM_ARGS][30];  /* command line arguments */
  38. char gnubuf[100];      /* command which is sent to gnuplot */
  39. char syscmd[100];      /* command which is sent to system */
  40. char filename[80];     /* name of data file */
  41. char topic[20];        /* help topic */
  42. int *co_lista;         /* lista for writing covariant matrix */
  43. double **alpha;        /* matrix and vectors used in LM algorithm */
  44. double *beta, *da;
  45. int num_indep = 1;     /* number of independent variables */
  46. int last_ma;           /* ma for last function */
  47. int last_num_indep;    /* num_indep for last function */
  48. int linflag = 0;       /* is function linear */
  49. int failed;            /* did function fail? */
  50.  
  51. /* each function has a comment, this holds the comment of the */
  52. /* current function.  If function is plotted using gnuplot, comment */
  53. /* should be f(x) = .... , where the comment defines the function */
  54. /* as it would on the gnuplot command line */
  55. char comment[160]; 
  56.  
  57. /* mode in which file is opened for reading and writing parameters */
  58. char mode[5]; 
  59. /* see description of data_order in fit.h */
  60. /* note that the order is the assignment of the data columns */
  61. /* internally. */
  62. struct data_order order;
  63. FILE *stream;         /* stream for reading and writing parameters to a file */
  64. int datarows = 1024;  /* number of rows in the data matrix */
  65. int datacols = 8;     /* number of columns in the data matrix */
  66. int echo = 0;
  67.  
  68. if(argc > 1) debug = atoi(argv[1]);
  69.  
  70. help("NOTE",100);
  71.  
  72. /* assign the proper command string to GNUPLOT */
  73.  
  74. strcpy(GNUPLOT,"gnuplot");
  75.  
  76. #ifdef OS2 
  77. strcpy(GNUPLOT,"gnuplot 2> gnuout");
  78. #endif
  79.  
  80. #ifdef DOS
  81. grflag = 0;
  82. #endif
  83.  
  84. /* under UNIX, this rather complicated command string is needed. */
  85. /* If modified, you should comment this one out and make another */
  86. /* so you can go back to this one if needed. */
  87. /* Things work pretty well with just a "gnuplot" command, */
  88. /* except that some of the things gnuplot writes to stderr */
  89. /* make it to the user's screen.  Simply using "gnuplot >& gnuout" */
  90. /* does not work, because popen() uses sh commands, and  */
  91. /* "gnuplot >& gnuout" uses csh redirection. */
  92.  
  93. #ifdef UNIX 
  94. strcpy(GNUPLOT,"/bin/csh -f -c \"gnuplot >& /tmp/gnuout\""); 
  95. /* strcpy(GNUPLOT,"gnuplot"); */
  96. if(debug)printf("*%s*\n", GNUPLOT);
  97. #endif
  98.  
  99. /* set up for one independent variable by default */
  100. order.x = ivector(num_indep);
  101. order.x[0] = 0;
  102. order.xsig = ivector(num_indep);
  103. for(i = 0; i < num_indep; i++) order.xsig[i] = -1;
  104.  
  105. xmin = dvector(num_indep);
  106. xmax = dvector(num_indep);
  107.  
  108. last_num_indep = num_indep;
  109.  
  110. /* assign these pointers to point to NULL */
  111. a=NULL;
  112. covar=NULL;
  113. lista=NULL;
  114. func=NULL;
  115. pstream = NULL;
  116.  
  117. last_ma = ma = 0;
  118.  
  119. /* allocate space for data matrix */
  120. data=dmatrix(datacols,datarows);
  121.  
  122. /* command processing loop */
  123.  
  124. while(strncmp(command,"quit",4) != 0){
  125.     strcpy(command,"\x0\x0\x0\x0\x0\x0\x0\x0");
  126.     printf("fit2> ");
  127.     gets(command);
  128.         /* if there is no command, reprompt */
  129.      if(strcmp(command,"")==0) continue; 
  130.  
  131.     if(echo) printf("%s\n",command);
  132.  
  133.         /* gd stands for get data */
  134.     if(strncmp(command,"gd",2) == 0){
  135.         if(debug) printf("in main, calling get_data()\n");
  136.         ndata=get_data(data,command,num_indep,inbuf,&order,filename,
  137.                             datarows, datacols);
  138.         if(debug) printf("in main, get_data() returned %d\n", ndata);
  139.         if(ndata ==0) printf("NO DATA, check filename\n");
  140.         printf("gd: %d data points\n", ndata);
  141.     }
  142.  
  143.         /* md stands for make data */
  144.     else if(strncmp(command,"md",2) == 0){
  145.         if(debug) printf("in main, calling make_data()\n");
  146.         failed=make_data(func, num_indep, a, ma, command,inbuf);
  147.         if(debug) printf("in main, make_data() returned %d\n", failed);
  148.     }
  149.  
  150.        /* sh stands for show */
  151.     else if(strncmp(command,"sh",2) == 0){
  152.         for(i = 0;  i < num_indep; i++)printf("order.x[%d]: %d\n", i, order.x[i]);
  153.         printf("order.y: %d\n", order.y);
  154.         printf("order.yfit: %d\n", order.yfit);
  155.         printf("order.sig: %d\n", order.sig);
  156.         printf("order.nsig: %d\n", order.nsig);
  157.         printf("order.ssig: %d\n", order.ssig);
  158.         printf("order.isig: %d\n", order.isig);
  159.         printf("order.osig: %d\n", order.osig);
  160.         for(i = 0;  i < num_indep; i++)printf("order.xsig[%d]: %d\n", i, order.xsig[i]);
  161.         if(func != NULL){
  162.             printf("function: %s %s\n", fname, comment);
  163.             printf("varying parameters: ");
  164.             for(i = 0; i < mfit; i++) printf(" %d",lista[i]);
  165.             printf("\n");
  166.         }
  167.         if(ndata) printf("data file: %s # data points: %d\n", filename, ndata);
  168.         if(wiflag){
  169.             printf("\n windowing on.");
  170.             for(i = 0; i < num_indep; i++) 
  171.                 printf("xmin%d: %g xmax%d: %g\n", i, xmin[i], i, xmax[i]);
  172.             printf("\n");
  173.         }
  174.         else printf("windowing off\n");
  175.         if(grflag)printf("graphing turned on\n");
  176.         else printf("graphing turned off\n");
  177.     }
  178.  
  179.         /* fn command tells us which function to fit to */
  180.     else if(strncmp(command,"fn",2) == 0){
  181.         last_ma = ma;
  182.         last_num_indep = num_indep;
  183.         iopt2 = sscanf(command,"%s %s %d ", inbuf, fname, &iopt1);
  184.         if(debug)printf("in main, calling getfcnptr()\n");
  185.         if((func = (void *)getfcnptr(fname,&num_indep, &linflag,&ma, comment)) != NULL){
  186.             if(debug)printf("in main, returned from getfcnptr()\n");
  187.  
  188.             /* get ma from command line or prompt for it if needed */
  189.             if( ma == -1){
  190.                 if(iopt2 > 2)
  191.                     ma = iopt1;
  192.                 else{
  193.                     printf("\n Enter number of parameters: ");
  194.                         gets(inbuf);
  195.                         sscanf(inbuf,"%d",&ma);
  196.                 }
  197.             }
  198.             if(debug)("ma: %d last_ma: %d", ma, last_ma);
  199.  
  200.             /* if covar, a, and lista are allocated and need to change size, we free them */
  201.             if(covar != NULL && ma != last_ma){
  202.                 if(debug)printf("free_dmatrix(covar,ma,ma)\n");
  203.                 free_dmatrix(covar,last_ma,last_ma);
  204.             }
  205.             if(lista !=NULL && ma != last_ma){
  206.                 if(debug)printf("free(lista)\n");
  207.                 free(lista);
  208.             }
  209.             if(a != NULL && ma != last_ma){
  210.                 if(debug)printf("free(a)\n");
  211.                 free(a);
  212.             }
  213.  
  214.             /* if num_indep changed, we reallocate some stuff */
  215.             if(num_indep != last_num_indep){
  216.                 free(order.x);
  217.                 free(order.xsig);
  218.                 free(xmin);
  219.                 free(xmax);
  220.                 order.x = ivector(num_indep);
  221.                 order.xsig = ivector(num_indep);
  222.                 xmin = dvector(num_indep);
  223.                 xmax = dvector(num_indep);
  224.  
  225.                 /* assign default values */
  226.                 for(i = 0; i < num_indep; i++){
  227.                     order.xsig[i] = -1;
  228.                     order.x[i] = i;
  229.                 }
  230.             }
  231.  
  232.             /* allocate space for covar, lista, and a */
  233.             if(ma != last_ma){
  234.                 covar=dmatrix(ma,ma);
  235.                 lista=ivector(ma);
  236.                 a=dvector(ma);
  237.                 /* initialize a's to 5, default to fitting all of the parameters */
  238.                 for(i=0;i<ma;i++){
  239.                     lista[i] = i;
  240.                     a[i] = 5;
  241.                 }
  242.                 mfit=ma;
  243.                 printf("%s: %d parameters\n", command, ma);
  244.             }
  245.         }
  246.         else printf("Function %s not found, try lf to list functions\n", fname);
  247.     }
  248.  
  249.         /* fit command does the fit */
  250.     else if(strncmp(command,"fit",2) == 0){
  251.         if(func == NULL)printf("fi failed, choose function first\n");
  252.         else if(ndata == 0)printf("fi failed, no data\n");
  253.         else{
  254.             sscanf(command,"%s %d ", inbuf, &itmax);
  255.             /* nonzero returns 1 if all sigma's are non-zero */
  256.             if(nonzero(data[order.sig], ndata)){
  257.                 if(debug)printf("in main, calling mrqfit(), itmax: %d\n", itmax);
  258.                 failed = mrqfit(data,order,num_indep,ndata,itmax,a,ma,
  259.                         lista,mfit,covar,&chisq,func,filename,comment);
  260.                 if(debug)printf("in main, mrqfit() returned %d\n", failed);
  261.             }
  262.             else{
  263.                 failed = 1;
  264.                 printf("all sigmay's must be non-zero, check weighting and order\n");
  265.             }
  266.             if(failed == 0)printf("fit done\n");
  267.             if(failed != 0)printf("fit failed\n");
  268.             printf("%s\n",comment);
  269.         }
  270.     }
  271.  
  272.         /* ip command initializes the parameters */
  273.     else if(strncmp(command,"ip",2) == 0){
  274.         if(debug)printf("in main, calling ip()\n");
  275.         ip(command, inbuf, func, a, ma);
  276.         if(debug)printf("in main, returned from ip()\n");
  277.     }
  278.  
  279.         /* cp command changes one of the parameters */
  280.     else if(strncmp(command,"cp",2) == 0){
  281.         if(debug)printf("in main, calling ip()\n");
  282.         cp(command, inbuf, func, a, ma);
  283.         if(debug)printf("in main, returned from ip()\n");
  284.     }
  285.  
  286.         /* sp command selects the parameters */
  287.     else if(strncmp(command,"sp",2) == 0){
  288.         if(debug)printf("in main, calling sp()\n");
  289.         mfit = sp(command, inbuf, func, lista, ma);
  290.         if(debug)printf("in main, sp() returned %d\n", mfit);
  291.     }
  292.  
  293.         /* pp command prints the parameters */
  294.     else if(strncmp(command,"pp",2) == 0){
  295.         if(func == NULL)
  296.             printf("pp failed, you must select a function first\n");
  297.         else{
  298.             printf("%s\n",comment);
  299.             for(i = 0; i < ma; i++) printf("a%d= %g\n", i, a[i]);
  300.             printf("chisqr = %g\n", chisq);
  301.             printf("%s\n",command);
  302.         }
  303.     }
  304.  
  305.         /* wp command writes the parameters to a file */
  306.     else if(strncmp(command,"wp",2) == 0){
  307.         if(func == NULL)
  308.             printf("wp failed, you must select a function first\n");
  309.         else{
  310.             strcpy(file,"a.dat");  /* filename defaults to a.dat */
  311.             if(sscanf(command,"%s %s %s", inbuf, file,mode) < 3) 
  312.                 strcpy(mode,"wt");  /* mode defaults to overwriting */
  313.             stream = fopen(file,mode);
  314.             if(stream){
  315.                 for(i=0; i < ma; i++)
  316.                     fprintf(stream,"%19.14e\n", a[i]);
  317.                 fclose(stream);
  318.             }
  319.             else
  320.                 printf("cannot open file %s in mode %s\n", file, mode);
  321.             printf("%s\n",command);
  322.         }
  323.     }
  324.  
  325.         /* wf writes the fitting function data to a file */
  326.     else if(strncmp(command,"wf",2) == 0){
  327.         if(func == NULL)
  328.             printf("wf failed, you must select a function first\n");
  329.         else{
  330.             if(debug)printf("in main, calling calc_yfit()\n");
  331.             failed = calc_yfit(func, data, order, num_indep, a, ndata, ma, &chisq);
  332.             if(debug)printf("in main, calc_yfit() returned %d\n", failed);
  333.             if(!failed){
  334.                 strcpy(file,"fit.dat");  /* filename defaults to fit.dat */
  335.                 sscanf(command,"%s %s", inbuf, file);
  336.                 stream = fopen(file,"wt");
  337.                 if(stream){
  338.                     for(i=0; i < ndata; i++){
  339.                         for(j = 0; j < num_indep; j++)
  340.                             fprintf(stream,"%g ", data[order.x[j]][i]);
  341.                         fprintf(stream,"%g\n", data[order.yfit][i]);
  342.                     }
  343.                     fclose(stream);
  344.                 }
  345.                 else
  346.                     printf("cannot open file %s\n", file);
  347.             }
  348.             else printf("%s failed \n",command);
  349.         }
  350.     }
  351.  
  352.         /* rp reads parameters from a file */
  353.     else if(strncmp(command,"rp",2) == 0){
  354.         if(func == NULL)
  355.             printf("rp failed, you must select a function first\n");
  356.         else{
  357.             strcpy(file,"a.dat");  /* filename defaults to a.dat */
  358.             sscanf(command,"%s %s", inbuf, file);
  359.             stream = fopen(file,"rt");
  360.             i = 0;
  361.             if(stream){
  362.                 while(fgets(inbuf,28,stream) && i < ma){
  363.                     a[i] = atof(inbuf);
  364.                     i++;
  365.                 }
  366.                 fclose(stream);
  367.                 printf("%s: %d parameters read in\n", command, i);
  368.             }
  369.             else printf("cannot open file %s\n", file);
  370.         }
  371.     }
  372.  
  373.         /* plot opens a pipe to gnuplot and tells it to plot the data and fit */
  374.     else if(strncmp(command,"plot",2) == 0){
  375.         if(func == NULL)
  376.             printf("plot failed, you must select a function first\n");
  377.         else{
  378.             grflag = 1;
  379.             if(debug)printf("in main, calling plot()\n");
  380.             failed=plot(func, data, order, num_indep, ndata, 
  381.                 filename, comment, a, ma);
  382.             if(debug)printf("in main, plot() returned %d\n", failed);
  383.             printf("%s\n",command);
  384.             if(failed)printf("%s failed\n",command);
  385.         }
  386.     }
  387.  
  388.         /* pr opens a pipe to gnuplot and tells it to plot the residual error */
  389.         /* "pr 1" plots bestfit vs y, "pr 2" plots (bestfit-y) vs x */
  390.     else if(strncmp(command,"pr",2) == 0){
  391.         if(func == NULL)
  392.             printf("pr failed, you must select a function first\n");
  393.         else{
  394.             grflag = 1;
  395.             if(debug)printf("in main, calling pr()\n");
  396.             failed=pr(command,func, data, order, num_indep, ndata, 
  397.                 filename, comment, a, ma);
  398.             if(debug)printf("in main, pr() returned %d\n", failed);
  399.             if(failed)
  400.                 printf("%s failed\n",command);
  401.         }
  402.     }
  403.  
  404.     else if(strncmp(command,"help",1) == 0){
  405.          strcpy(topic,"");
  406.         sscanf(command,"%*s %s", topic);
  407.         if(debug)printf("in main, calling help()\n");
  408.         help(topic,100);
  409.         if(debug)printf("in main, returned from help()\n");
  410.     }
  411.  
  412.         /* send commands beginning with set or load to gnuplot */
  413.     else if(strncmp(command,"set",2) == 0 ||  strncmp(command,"lo",2) == 0){
  414.     sprintf(gnubuf,"%s\n",command);
  415.     gnucmd(gnubuf);
  416.     }
  417.  
  418.         /* wt select weighting for fit */
  419.     else if(strncmp(command,"wt",2) == 0){
  420.         if(ndata == 0)
  421.             printf("wt failed, you must get data first\n");
  422.         else{
  423.             sscanf(command,"%*s %s", inbuf);
  424.             if(strncmp(inbuf,"statistical",1) == 0){
  425.                 order.sig = order.ssig;
  426.                 for(i=0; i<ndata; i++)
  427.                     data[order.sig][i] = sqrt(fabs(data[order.y][i]));
  428.             }
  429.             else if(strncmp(inbuf,"instrumental",1) == 0)
  430.                 order.sig = order.isig;
  431.             else if(strncmp(inbuf,"none",1) == 0)
  432.                 order.sig = order.nsig;
  433.             else if(strncmp(inbuf,"other",1) == 0)
  434.                 order.sig = order.y;
  435.             else help("wt",1);
  436.         }
  437.     }
  438.  
  439.         /* or redefines order of colums in data */
  440.     else if(strncmp(command,"or",2) == 0){
  441.         if(ndata == 0)
  442.             printf("or failed, you must get data first\n");
  443.         else{
  444.             i = parse(command, cmdargs);
  445.             if(debug) printf("num_indep: %d num args: %d\n", num_indep, i);
  446.             for(j = 0; j < num_indep; j++) if(i > j){
  447.                 if(debug) printf("cmdargs[%d]: %s\n",j, cmdargs[j]);
  448.                 order.x[j] = atoi(cmdargs[j]);
  449.                 if(debug) printf("order.x[%d]: %d\n", j, order.x[j]);
  450.             }
  451.             if(i > num_indep) order.y = atoi(cmdargs[num_indep]);
  452.             if(i > num_indep + 1) order.isig = atoi(cmdargs[num_indep + 1]);
  453.             for(j = 0; j < num_indep; j++)
  454.                 if(i > j+num_indep+2) order.xsig[j] = atoi(cmdargs[j+num_indep+2]);
  455.             if(i == 0) help("or",1);
  456.         }
  457.     }
  458.  
  459.         /* co prints the covariant matrix from the fit */
  460.         /* if there are no arguments, output is to */
  461.         /* screen.  An argument is taken as a filename */
  462.         /* and matrix is written to a file */
  463.     else if(strncmp(command,"co",2) == 0){
  464.         if(func == NULL)
  465.             printf("co failed, you must select a function first\n");
  466.         else{
  467.             strcpy(file,"");
  468.             sscanf(command,"%s %s", inbuf, file);
  469.             co_lista = ivector(ma);
  470.             beta = dvector(ma);
  471.             da = dvector(ma);
  472.             alpha=dmatrix(ma,ma);
  473.             for(i = 0; i < ma;  i++) co_lista[i] = i;
  474.             if(debug)printf("in main, calling alpha_beta_chisqr()\n");
  475.             alpha_beta_chisq(data, order, num_indep, ndata,a,ma,
  476.                                     co_lista,ma,alpha,beta,&chisq,func);
  477.             if(debug)printf("in main, returned from alpha_beta_chisqr()\n");
  478.             if(debug)printf("in main, calling solve_for_da()\n");
  479.             solve_for_da(alpha, covar, beta, da, ma);
  480.             if(debug)printf("in main, returned from solve_for_da()\n");
  481.             free(co_lista);
  482.             free(beta);
  483.             free(da);
  484.             free_dmatrix(alpha,ma,ma);
  485.             if( strcmp(file,"") == 0){
  486.                 for(i = 0; i < ma; i++){
  487.                     strcpy(inbuf,"");
  488.                     for(j = 0; j < ma; j++){
  489.                         sprintf(buf,"%9.2g ",covar[i][j]);
  490.                         strcat(inbuf, buf);
  491.                     }
  492.                     printf("%s\n",inbuf);
  493.                 }
  494.             }
  495.             else{
  496.                 stream = fopen(file,"wt");
  497.                 for(i = 0; i < ma; i++){
  498.                     strcpy(inbuf,"");
  499.                     for(j = 0; j < ma; j++){
  500.                         sprintf(buf,"%10g ",covar[i][j]);
  501.                         strcat(inbuf, buf);
  502.                     }
  503.                 fprintf(stream,"%s\n",inbuf);
  504.                 }
  505.                 fclose(stream);
  506.             }
  507.         }
  508.     }
  509.  
  510.         /* ad allocates data matrix to a specified size */
  511.     else if(strncmp(command,"ad",2) == 0){
  512.         iopt1 = 0;
  513.         iopt2 = 0;
  514.         sscanf(command,"%s %d %d", inbuf, &iopt1, &iopt2);
  515.         if(debug)printf("cols: %d rows: %d\n", iopt1, iopt2);
  516.         if(iopt1 > 3 && iopt2 > 1 && iopt1 < 128){
  517.             if(debug) printf("calling free_dmatrix(data,datacols,datarows);\n");
  518.             free_dmatrix(data,datacols,datarows);
  519.             datacols = iopt1;
  520.             datarows = iopt2;
  521.             if(debug) printf("calling data=dmatrix(data,datacols,datarows);\n");
  522.             data=dmatrix(datacols,datarows);
  523.             if(debug) printf("returned from data=dmatrix(data,datacols,datarows);\n");
  524.             ndata = 0;
  525.         }
  526.         else printf("ad: columns rows, columns => columns in data file + 3\n");
  527.     }
  528.  
  529.         /* gr turns graphing on and off */
  530.     else if(strncmp(command,"gr",2) == 0){
  531.         sscanf(command,"%*s %d", &iopt1);
  532.         if(iopt1 == 1|| iopt1 == 0){
  533.             grflag =iopt1;
  534.         }
  535.         else printf("must be \"gr 0\" (graphing off) or \"gr 1\" (graphing on)\n");
  536.     }
  537.  
  538.         /* lf lists the functions available for fitting */
  539.     else if(strncmp(command,"lf",2) == 0){
  540.         listfcns();
  541.     }
  542.  
  543.         /* wi with one parameter turns windowing on and off */
  544.         /* with two parameters, it selects xmin and xmax */
  545.     else if(strncmp(command,"wi",2) == 0){
  546.         iopt1 = parse(command,cmdargs);
  547.         if(iopt1 < 1) printf("wi failed, try help wi\n");
  548.         else if(iopt1 == 1){
  549.             wiflag = atoi(cmdargs[0]);
  550.         }
  551.         else if(iopt1 > 2*num_indep) printf("wi failed, try help wi\n");
  552.         else{
  553.             if(debug)printf("%d arguments\n", iopt1);
  554.             wiflag = 1;
  555.             i = 0;
  556.             while(i <= iopt1 && (i/2) < num_indep){
  557.                 xmin[i/2] = atof(cmdargs[i]);
  558.                 xmax[i/2] = atof(cmdargs[i + 1]);
  559.                 i +=2;
  560.             }
  561.         }
  562.     }
  563.  
  564.         /* ve selects verbosity */
  565.     else if(strncmp(command,"ve",2) == 0){
  566.         if(sscanf(command,"%*s %d",&iopt1) > 0 && ( iopt1 ==2 || iopt1 == 1|| iopt1 == 0)){
  567.             veflag =iopt1;
  568.         }
  569.         else printf("must be ve 0 1 or 2\n");
  570.     }
  571.  
  572.         /* pa pauses for a number of seconds */
  573.     else if(strncmp(command,"pause",2) == 0){
  574.         if(sscanf(command,"%*s %d",&iopt1) > 0 && ( iopt1 > 0)){
  575.             sleep(iopt1);
  576.         }
  577.         else help("pause",1);
  578.     }
  579.  
  580.        /* li does a linear least squares fit */
  581.     else if(strncmp(command,"li",2) == 0){
  582.     if(func == NULL){
  583.         printf("You must select a function first\n");
  584.         failed = 1;
  585.         }
  586.     else if(ndata < 1){
  587.         printf("No data\n");
  588.         failed = 1;
  589.     }
  590.     else if(linflag ==0){
  591.         printf("Not a linear function\n");
  592.         failed = 1;
  593.     }
  594.     else if(nonzero(data[order.sig], ndata)){
  595.         if(debug)printf("in main, calling linear_fit()\n");
  596.         failed = linear_fit(data,order,num_indep,ndata,itmax,a,ma,
  597.                         lista,mfit,covar,&chisq,func,filename,comment);
  598.         if(debug)printf("in main, linear_fit() returned %d\n", failed);
  599.     }
  600.     else{
  601.         failed = 1;
  602.         printf("all sigmay's must be non-zero, check weighting and order\n");
  603.     }
  604.     if(failed == 0){
  605.         printf("fit done\n");
  606.         printf("%s\n",comment);
  607.     }
  608.     if(failed != 0)printf("fit failed\n");
  609.     }
  610.  
  611.         /* de selects debugging option */
  612.     else if(strncmp(command,"de",2) == 0){
  613.         if(sscanf(command,"%*s %d",&iopt1) > 0 && ( iopt1 ==2 || iopt1 == 1|| iopt1 == 0)){
  614.             debug =iopt1;
  615.         }
  616.         else printf("must be de 0 1 or 2\n");
  617.     }
  618.  
  619.         /* echo selects debugging option */
  620.     else if(strncmp(command,"echo",2) == 0){
  621.         if(sscanf(command,"%*s %d",&iopt1) > 0 && (iopt1 == 1|| iopt1 == 0)){
  622.             echo =iopt1;
  623.         }
  624.         else printf("must be echo 0 or 1\n");
  625.     }
  626.  
  627.     else if(strncmp(command,"run",2) == 0){
  628.         if((i = parse(command, cmdargs)) > 0){
  629.             strcpy(syscmd,"");
  630.             for(j = 0; j < i; j++){
  631.                 strcat(syscmd, cmdargs[j]);
  632.                 strcat(syscmd," ");
  633.                 if(debug) printf("j: %d syscmd *%s* \n",j,syscmd);
  634.             }
  635.             if(debug) printf("Sending command %s to system\n",syscmd);
  636.             system(syscmd);
  637.         }
  638.     }
  639.  
  640.         /* gn send commands to gnuplot */
  641.     else if(strncmp(command,"gn",2) == 0){
  642.         if((i = parse(command, cmdargs)) > 0){
  643.             strcpy(gnubuf,"");
  644.             for(j = 0; j < i; j++){
  645.                 strcat(gnubuf, cmdargs[j]);
  646.                 strcat(gnubuf," ");
  647.                 if(debug) printf("j: %d gnubuf *%s* \n",j,gnubuf);
  648.             }
  649.             sprintf(gnubuf,"%s\n",gnubuf);
  650.             gnucmd( gnubuf);
  651.         }
  652.         else help("gn",1);
  653.     }
  654.  
  655.     else if(strncmp(command,"quit",4))
  656.         printf("command %s not recognized\n",command);
  657. }
  658.  
  659.  
  660. /* free allocated arrays and close open pipe to gnuplot */
  661. free_dmatrix(data,datacols,datarows);
  662. if(a != NULL) free(a);
  663. if(covar != NULL) free_dmatrix(covar,ma,ma);
  664. if(lista != NULL) free(lista);
  665. free(order.x);
  666. free(order.xsig);
  667. free(xmin);
  668. free(xmax);
  669. #ifndef DOS
  670. if(pstream != NULL){
  671.     fprintf(pstream,"quit\n");
  672.     pclose(pstream);
  673. }
  674. #endif
  675. #ifdef UNIX
  676. if(!debug){
  677.     system("rm /tmp/fit.tmp");
  678.     system("rm /tmp/gnuout");
  679. }
  680. #endif
  681. printf("\n");
  682. return 0;
  683. }
  684.  
  685. int gnucmd(char command[100]){
  686.  
  687. #ifndef DOS
  688.  
  689. grflag = 1;
  690. /* open the pipe */
  691. if(pstream == NULL){
  692.     if(debug) printf("attempting to open pipe to gnuplot\n");
  693.     pstream = popen(GNUPLOT,"w");  
  694. }    
  695. if(pstream == NULL){
  696.     printf("Failed to open pipe to gnuplot\n");
  697.     return 1;
  698. }
  699. fprintf(pstream, command);
  700. if(debug) printf("command to gnuplot: %s", command);
  701. fflush(pstream);
  702. #endif
  703.  
  704. #ifdef DOS
  705. printf("The DOS version of this program does not support gnuplot\n");
  706. printf("directly.  You need a real operating system for that.  \n");
  707. printf("I suggest OS/2 2.1. You can use wf to write the fit  \n");
  708. printf("to a file and use gnuplot by itself for plotting  \n\n");
  709. grflag = 0;
  710. #endif
  711.  
  712. return 0;
  713. }
  714.