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