#include "xstcl.h"
#include "xscf.h"
#include <math.h>



#define MAXOPTS 41




int
xs_tclout(ClientData cdata,Tcl_Interp* interp,int objc, Tcl_Obj* CONST objv[] )
{
        static char* outToTcl;
        char* chosen = NULL;
        char* next = NULL;
/*
        pointers for xspec common block variables

*/
        float* energies = NULL;         /* ear */
        float* phaCor = NULL;   
        float* phaObs = NULL;
        float* phaMod = NULL; 
        float* phaBkg = NULL;
        float* varObs = NULL;
        float* varBkg = NULL;
        float* areasc = NULL;
	float* backsc = NULL;
        float* effArea = NULL;          /* effar  */
        float* flux  = NULL;            /* fluxar */
        float* components = NULL;       /* esavar */
        float* responseElts = NULL;     /* respon */
        int* group = NULL;              /* igroup */
        int* beginChan = NULL;          /* ichanb */
        int* endChan = NULL;            /* ichane */
        double* values = NULL;

        int logging = 0;
        int status  = 0;
        int retStatus  = TCL_OK;
        int posn = 0;
        int plot = 0;
        int interactive = 0;
        char* tmp=NULL;
        int paramNo = 1;
        int componentNo = 1;
        int dataNo = 1;
        int nparams;
        int i = 0;
	int source = 1;
	float energy = 0.;
        float deltaEnergy = 0.;
        char* tmp1 = NULL;
        char* tmp2 = NULL;        
        unsigned int paramRequested = 1; /* for genfit */

        unsigned int selected = 1;
        
        const char* option[] = {
                "",             /* [0]  null selection - put there to make 0 an invalid input
                                         condition */
		"areascal",     /* [1] areascal values */
		"arf",          /* [2] arf filename */
		"backgrnd",     /* [3] background filename */
		"backscal",     /* [4] backscal values */
                "chatter",      /* [5]  current chatter level */
                "compinfo",     /* [6]  compinfo <comp.no.>   -> the name, start model param # 
                                         and number of parameters of the specified model 
                                         component */
		"cosmo",        /* [7]  cosmological parameters in use */
                "datagrp",      /* [8]  number of datagroups */
                "datasets",     /* [9]  number of datagroups */
                "dof",          /* [10]  no. of d.o.f. */
                "energies",     /* [11]  energies <dataset no.> -> response energies for specified 
                                         dataset */
                "eqwidth",      /* [12]  eqwidth <datagroup no.> -> last equivalent width 
                                         calculated for specified datagroup */
                "error",        /* [13]  error <par.no.>  -> last confidence region calculated 
                                        for specified parameter */
		"expos",        /* [14]  exposure time */
                "filename",     /* [15]  filename <dataset no.> -> filename of specified dataset */
                "flux",        /* [16]  flux <dataset no.> -> last model (energy flux, low range on error, high range on error)
                                         calculated for specified dataset */
		"ftest",        /* [17] last calculated F statistic */
                "genpop",       /* [18]  genpop <par.no.> -> the genetic algorithm population for 
                                         specified parameter */
		"goodness",     /* [19] last calculated goodness */
                "idline",       /* [20] idline <energy> <delta> -> get possible line IDs in 
                                         specified energy range */
                "lumin",        /* [21]  luminosity <dataset no.> -> last model (luminosity, low range on error, high range on error)
                                         calculated for specified dataset */
                "model",       /* [22] description of current model */
                "modcomp",      /* [23] number of model components */
                "modpar",       /* [24] number of model parameters */
                "modval",       /* [25] modval   <dataset no.>  current calculated model values 
                                         for specified dataset */
                "noticed",      /* [26] noticed <dataset no.>  (low,high) noticed channels for 
                                         specified dataset */
                "param",        /* [27] param <par.no.> -> (initial,delta,min,low,high,max) for
                                         specified model parameter */
                "peakrsid",     /* [28] peakrsid <dataset no.> -> energies and strengths of the 
                                         peak residuals (+ve and -ve) for the specified dataset. 
                                         If an additional two arguments are given, they  are an 
                                         energy range in which to search */
                "pinfo",        /* [29] pinfo <par.no.> -> parameter name and unit */
                "plink",        /* [30] plink <par.no.> -> information on parameter linking. 
                                         (T/F, parameter, multiplicative factor, 
                                         additive constant)*/
                "plot",        /* [31] plot <option> <array> <plot group no.> -> the output from 
                                         a plot command. The array options are x/xerr/y/yerr */
                "plotgrp",      /* [32] number of plot groups */
		"query",        /* [33] state of the query switch */
                "rate ",        /* [34] rate <dataset no.> -> count rate for specified dataset */
		"response",     /* [35] response filename */
		"simpars",      /* [36] simulated model parameters */
                "solab",        /* [37] current solar abundance table values */
                "stat",         /* [38] value of statistic */
                "varpar",       /* [39] number of variable fit parameters */
		"xflt",         /* [40] value of XFLT#### keywords */
                "?"             /* [41] list of valid options */
        };
              
        /* minimum number of characters in the string required to identify uniquely */
        const int optLen[] = {0, 3, 3, 5, 5, 2, 3, 3, 5, 5, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
        /*                    0  1  2  3  4  5  6  7  8  9  1  1  2  3  4  5  6  7  8  9*/
                              1, 1, 4, 4, 4, 4, 1, 2, 2, 2, 3, 5, 5, 1, 2, 2, 2, 2, 2, 1, 1, 1} ; 
        /*                    2  1  2  3  4  5  6  7  8  9  3  1  2  3  4  5  6  7  8  9  4  1*/    
       const char* celt[] = {"H ", "He", "C ", "N ", "O ", "Ne", "Na", "Mg","Al", "Si", "S ",
                 "Cl", "Ar", "Ca", "Cr", "Fe","Co", "Ni"};

        char* inString = NULL;
        Tcl_Obj* inputObj = NULL;
        int objlen = 0;
        int arglen = 0;
        inString = (char *)malloc(256*sizeof(char));
        
        inputObj = Tcl_ConcatObj(objc-1,objv+1);
        
        strcpy (inString, Tcl_GetStringFromObj(inputObj,&objlen));

        
        
 
	if ( objv[1] == NULL || objlen == 0)
	{
		selected = 0;
	}
	else
	{
 		int j = 0;
                char* retrieved = NULL;
		retrieved = Tcl_GetStringFromObj(objv[1],&j);
                chosen = (char *)calloc(j+1,sizeof(char));
		
                strcpy(chosen,retrieved);

        	while (selected <= (unsigned int)MAXOPTS )
                {
                        if (!STRNCASECMP(chosen,option[selected],optLen[selected])) break;
                        ++selected;
                }
        }
        
        if (selected == 0 || selected >= MAXOPTS)
        {
                /* didn't find the selection */
                char* errStrg= NULL;
                int iout = 1;
                char* allopts = NULL;
                allopts = (char*)malloc((MAXOPTS*12 + 1)*sizeof(char));
                errStrg = (char*)malloc(72*sizeof(char));
                strcpy(allopts,"");
                if (selected == 0 ) 
                   sprintf(errStrg,"Syntax: TCLOUT <option> [args] - valid options are:\n");
                else 
                {
                        if (selected > MAXOPTS)  
                        {   
                                sprintf(errStrg,
                                "%s is not currently a valid TCLOUT option - valid options are\n",
                                chosen);
                        }
			else strcpy(errStrg,"Valid TCLOUT options are:\n");
                }
                for (iout = 1; iout <= MAXOPTS; ++iout)
                {
                        if (iout % 8 != 1 ) strcat(allopts,"\t");
                        strcat(allopts,option[iout]);
                        if ( iout % 8 == 0 ) strcat(allopts,"\n");
                }
                errStrg = realloc(errStrg,(strlen(errStrg)+strlen(allopts) +1)*sizeof(char));
                strcat(errStrg,allopts);
                XWRITE(errStrg,10);
                free(errStrg);
                free(allopts);
                if (selected > MAXOPTS) retStatus = TCL_ERROR;
		free(inString);
		return retStatus;
        }
        
        /* allocate for the output. The intention is to reallocate this string
           dynamically depending on requirements. A lot of the outputs are for an
           integer of width 8, so make this the default (remembering trailing null).*/
           
        outToTcl = (char*)realloc(outToTcl,9*sizeof(char));
        
       
/* we now have unsigned int 'selected' for use in a switch/case construction. The only complication is
   that some options take arguments that need to be processed further, and these are
   bunched together rather than treated in numerical order */

        switch (selected)
        {

                /* options that require dataset number */
         	case  1: /* areascal */
                case  2: /* arf */
                case  3: /* backgrnd */
                case  4: /* backscal */
                case 11: /* energies */
          	case 12: /* eqwidth */
          	case 14: /* expos */
                case 15: /* filename */
                case 16: /* flux */
	        case 21: /* luminosity */
                case 25: /* modval */
                case 26: /* noticed */
                case 28: /* peakrsid */
                case 34: /* rate */
                case 35: /* response */
                case 40: /* xflt */
                        if (numberOfDataSets() > 1)
                        {
                           if ( objc < 3)
                           {
                              char* errMessage = NULL;
                              errMessage = (char *)malloc(64*sizeof(char));
                              sprintf(errMessage,
                              "Dataset number argument required for option %s",chosen);
                              XWRITE(errMessage,10);
                              free(errMessage);
                              retStatus = TCL_ERROR;
                              break;
                                         
                           }
                           next = Tcl_GetStringFromObj(objv[2],&arglen);
                           if (!isInteger(next))
                           {
                              char errString[] = 
                                "Argument for dataset number is not an integer";
                              XWRITE(errString, 10);
                              retStatus = TCL_ERROR;
                              break;
                           }
                           sscanf(next,"%d",&dataNo);
                           if ( dataNo <= 0 || dataNo > numberOfDataSets())
                           {
                              char* errMessage = NULL;
                              errMessage = (char *)malloc(64*sizeof(char));
                              sprintf(errMessage,
                                "Invalid dataset number d: must satisfy 1 < d < %d ",
                                numberOfDataSets() );
                              XWRITE(errMessage,10);
                              free(errMessage);
                              retStatus = TCL_ERROR;
                              break;
                           }
                        }

/* for areascal, backscal, and expos we look for a source/back argument */

                        if (selected == 1 || selected == 4 || selected == 14 )
		        {
		      	   source = 1;
			   if ( objc >= 3 )
			   {
			      next = Tcl_GetStringFromObj(objv[2],&arglen);
			      if (!STRNCASECMP(next,"b",1)) source = 0;
                           }
			   if ( objc == 4 )
			   {
			      next = Tcl_GetStringFromObj(objv[3],&arglen);
			      if (!STRNCASECMP(next,"b",1)) source = 0;
                           }
			   if ( source == 0 && dataBkgTime(dataNo) == 0.0 )
			   {
                              char* errMessage = NULL;
                              errMessage = (char *)malloc(64*sizeof(char));
                              sprintf(errMessage,
				      "There is no background file for dataset %d",
                                      dataNo);
                              XWRITE(errMessage,10);
                              free(errMessage);
                              retStatus = TCL_ERROR;
                              return retStatus;
			   }
			}

/* now go through the options for this sub-case */

                        if (selected == 1) /* areascal */
			{
                                
			  if (source) {
			    areasc = udmArray("are");
			  } else {
			    areasc = udmArray("bre");
			  }

                          tmp = (char *)malloc((14+1)*sizeof(char));
                          outToTcl = realloc(outToTcl,((endBin(dataNo)-startBin(dataNo)+1)*14+1)*sizeof(char));

                          strcpy(outToTcl,"");

                          for (i = startBin(dataNo); i <= endBin(dataNo); i++)
			  {
                             sprintf(tmp," %13.5g",areasc[i-1]);
                             strcat(outToTcl,tmp);
                          }
                          free(tmp);
			}

                        if (selected == 2) /* arf */
                        {
                                outToTcl = realloc(outToTcl,256*sizeof(char));
                                sprintf(outToTcl,"%s",auxFilename(dataNo,"a"));
                                outToTcl = realloc(outToTcl,(strlen(outToTcl)+1)*sizeof(char));
                        }

                        if (selected == 3) /* background */
                        {
                                outToTcl = realloc(outToTcl,256*sizeof(char));
                                sprintf(outToTcl,"%s",auxFilename(dataNo,"b"));
                                outToTcl = realloc(outToTcl,(strlen(outToTcl)+1)*sizeof(char));
                        }

                        if (selected == 4) /* backscal */
			{
			  if (source) {
			    backsc = udmArray("reg");
			  } else {
			    backsc = udmArray("brg");
			  }
                                
                          tmp = (char *)malloc((14+1)*sizeof(char));
                          outToTcl = realloc(outToTcl,((endBin(dataNo)-startBin(dataNo)+1)*14+1)*sizeof(char));
                          strcpy(outToTcl,"");

                          for (i = startBin(dataNo); i <= endBin(dataNo); i++)
			  {
                             sprintf(tmp," %13.5g",backsc[i-1]);
                             strcat(outToTcl,tmp);
                          }
                          free(tmp);
			}

                        if (selected == 11) /* energies */
                        {
                                int lenOut = 1;
                                tmp = (char *)malloc((14+1)*sizeof(char));
                                strcpy(outToTcl,"");
                                
                                energies = udmArray("e");
                                
                                for ( i = firstEnergy(dataNo); i <= lastEnergy(dataNo); i++)
                                {
                                        sprintf(tmp," %13.5g",energies[i]);
                                        lenOut += 14;
                                        outToTcl = realloc(outToTcl,(lenOut+1)*sizeof(char));
                                        strcat(outToTcl,tmp);
                                }
                                free(tmp);
                        }

                        if (selected == 12) /* eqwidth */
                        {
                                outToTcl = realloc(outToTcl,(14+1)*sizeof(char));
                                sprintf(outToTcl," %13.5g",dataEqWidth(dataNo));
                        }

                        if (selected == 14) /* exposure */
			{
			    if (source) {
                                outToTcl = realloc(outToTcl,(14+1)*sizeof(char));
                                sprintf(outToTcl," %13.5g",dataSrcTime(dataNo));
			    } else {
                                outToTcl = realloc(outToTcl,(14+1)*sizeof(char));
                                sprintf(outToTcl," %13.5g",dataBkgTime(dataNo));
			    }
			}

                        if (selected == 15) /* filename */
                        {
                                outToTcl = realloc(outToTcl,256*sizeof(char));
                                sprintf(outToTcl,"%s",dataFileName(dataNo));
                                outToTcl = realloc(outToTcl,(strlen(outToTcl)+1)*sizeof(char));
                        }

                        if (selected == 16) /* flux */
                        {
                                flux  = udmArray("f");
                                outToTcl = realloc(outToTcl,(42+1)*sizeof(char));
                                sprintf(outToTcl," %13.5g %13.5g %13.5g",dataFlux(dataNo),dataFluxErr(dataNo,"l"),dataFluxErr(dataNo,"h"));
                        }

                        if (selected == 21) /* luminosity */
                        {
                                flux  = udmArray("f");
                                outToTcl = realloc(outToTcl,(42+1)*sizeof(char));
                                sprintf(outToTcl," %13.5g %13.5g %13.5g",dataLumin(dataNo),dataLuminErr(dataNo,"l"),dataLuminErr(dataNo,"h"));
                        }

                        if (selected == 25) /* modval */
                        {
                                int lenOut = 1;
				int ii = 0;
                                tmp = (char *)malloc((14+1)*sizeof(char));
                                strcpy(outToTcl,"");
                                
                                for (ii = endPrevResp(dataNo); 
                                     ii <= endPrevResp(dataNo) + numberOfEnergies(dataNo); ii++)
                                {
                                	flux  = udmArray("f");
                                        sprintf(tmp," %13.5g",flux[ii]);
                                        lenOut += 14;
                                        outToTcl = realloc(outToTcl,(lenOut+1)*sizeof(char));
                                        strcat(outToTcl,tmp);
                                }
                                free(tmp);
                        }
                        if (selected == 26) /* noticed */
                        {
                                unsigned int inNoticedRange = 0;
                                unsigned int firstChan = 0;
                                unsigned int lastChan = 0;
                                int lenOut = 1;
                                tmp = malloc((14+1)*sizeof(char));
                                strcpy(outToTcl,"");
                                for (i=1; i <= dataGroupedChan(dataNo); ++i)
                                {
                                        if ( chanIsNoticed(dataNo,i) )
                                        {
                                                if (!inNoticedRange)
                                                {
                                                        firstChan = i;
                                                        inNoticedRange = 1;
                                                }
                                        }
                                        else
                                        {
                                                if (inNoticedRange)
                                                {
                                                        lastChan = i - 1;
                                                        sprintf(tmp," %5d - %5d",firstChan,
                                                          lastChan);
                                                        lenOut += 14;
                                                        outToTcl = 
                                                                realloc(outToTcl,(lenOut+1)*sizeof(char));
                                                        strcat(outToTcl,tmp);
                                                        inNoticedRange = 0;
                                                }
                                        }
                                }
                                if ( inNoticedRange )
                                {
                                        lastChan = dataGroupedChan(dataNo);
                                        sprintf(tmp," %5d - %5d",firstChan,lastChan);
                                        lenOut +=14;
                                        outToTcl = realloc(outToTcl,(lenOut+1)*sizeof(char));
                                        strcat(outToTcl,tmp);      
                                }
                                free(tmp);
                        }
                        if (selected == 28) /* peakrsid */
                        {
                                float peakE[2];
                                float peakNM[2];
                                /* this somewhat obscure setting tries to duplicate
                                   xparse's parse position as required in the fortran string.
                                   first, find the first blank character in the string. Then
                                   move to the next blank, using the strchr to locate the 
                                   first 'p' in the string (there must be one since if we
                                   are here the user typed tclout pe[akrsid]. The +1 is
                                   because C starts counting at 0. */
                                posn = strcspn(inString," ");
                                posn += strcspn(strchr(inString,'p')," ") + 1;
                                phaCor = udmArray("cor");
                                phaObs = udmArray("src");
                                phaBkg = udmArray("bkg");
                                phaMod = udmArray("mod");
                                energies = udmArray("e");
                                effArea  = udmArray("a");
                                peakResidual(inString,&posn,phaCor,phaObs,phaMod,phaBkg,
                                        energies,effArea,dataNo,peakE,peakNM);
                                outToTcl = realloc(outToTcl,(56+1)*sizeof(char));
                                sprintf(outToTcl," %13.5g %13.5g %13.5g %13.5g",
                                        peakE[0],peakNM[0],peakE[1],peakNM[1]);
                        
                        }
                        if (selected == 34) /* rate */
                        {
                                double accumA = 0.;
                                double accumB = 0.;
                                double accumC = 0.;
                                float corNorm = 0.;
                                phaObs = udmArray("src");
                                phaBkg = udmArray("bkg");
                                varObs = udmArray("osv");
                                varBkg = udmArray("obv");
                                areasc = udmArray("are");
                                phaMod = udmArray("mod");
                                for (i = startBin(dataNo); i <= endBin(dataNo); i++)
                                {
                                        accumA += phaObs[i-1]*areasc[i-1];
                                        accumB += varObs[i-1]*areasc[i-1]*areasc[i-1];
                                }
                                if ( dataBkgTime(dataNo) != 0)
                                {
                                        for (i = startBin(dataNo); i <= endBin(dataNo); i++)
                                        {
                                                accumA -= phaBkg[i-1]*areasc[i-1];
                                                accumB += varBkg[i-1]*areasc[i-1]*areasc[i-1];
                                        }
                                }
                                
                                if ( correctionIsActive(dataNo) )
                                {
                                	phaCor = udmArray("cor");
                                        corNorm = -dataCorNorm(dataNo);
                                        for (i = startBin(dataNo); i <= endBin(dataNo); i++)
                                        {
                                                accumA += corNorm*phaCor[i-1]*areasc[i-1];
                                        }
                                }
                                if ( numberOfComponents() > 0 )
				{
                                        for (i = startBin(dataNo); i <= endBin(dataNo); i++)
                                        {
                                                accumC += phaMod[i-1]*areasc[i-1];
                                        }
                                }

                                outToTcl = realloc(outToTcl,(45+1)*sizeof(char));
                                sprintf(outToTcl," %14.4g %14.4g %14.4g",accumA,
                                        sqrt(accumB), accumC);
                        }

                        if (selected == 35) /* response */
                        {
                                outToTcl = realloc(outToTcl,256*sizeof(char));
                                sprintf(outToTcl,"%s",auxFilename(dataNo,"r"));
                                outToTcl = realloc(outToTcl,(strlen(outToTcl)+1)*sizeof(char));
                        }

                        if (selected == 40) /* xflt */
			{
                           int lenOut = 4;
			   int ii = 0;
                           tmp = (char *)malloc(80*sizeof(char));
			   sprintf(tmp, " %3d", NumXFLT(dataNo));
                           strcpy(outToTcl,tmp);
                                
                           for (ii = 1; ii <= NumXFLT(dataNo); ii++)
                           {
                              sprintf(tmp," %s",getXFLT(ii,dataNo));
                              lenOut += strlen(tmp)+1;
                              outToTcl = realloc(outToTcl,lenOut*sizeof(char));
                              strcat(outToTcl,tmp);
                           }
                           free(tmp);

                        }

                        break;
                case  5: /* chatter */
                        XparseGetChatterLevel(&interactive,&logging); 
                        sprintf(outToTcl," %3d %3d", interactive, logging);
                        break;
                case  6: /* compinfo */
                        if (numberOfComponents() > 1)
                        {
                                if ( objc < 3)
                                {
                                        char* errMessage = NULL;
                                        errMessage = (char *)malloc(64*sizeof(char));
                                        sprintf(errMessage,
                                         "Component number argument required for option %s",chosen);
                                        XWRITE(errMessage,10);
                                        free(errMessage);
                                        retStatus = TCL_ERROR;
                                        break;
                                         
                                }

                                next = Tcl_GetStringFromObj(objv[2],&arglen);
                                if (!isInteger(next))
                                {
                                        char errString[] = 
                                           "Argument for component number is not an integer";
                                        XWRITE(errString, 10);
                                        retStatus = TCL_ERROR;
                                        break;
                                }
                                sscanf(next,"%d",&componentNo);
                                if ( componentNo <= 0 || componentNo > numberOfComponents())
                                {
                                        char* errMessage = NULL;
                                        errMessage = (char *)malloc(64*sizeof(char));
                                        sprintf(errMessage,
                                          "Invalid component number d: must satisfy 1 < d < %d ",
                                          numberOfComponents());
                                        XWRITE(errMessage,10);                                        
                                        free(errMessage);
                                        retStatus = TCL_ERROR;
                                        break;
                                }
                        }     
                        outToTcl = realloc(outToTcl,(24+1)*sizeof(char));
                        sprintf(outToTcl,"%8s  %6d  %6d",componentName(componentNo),
                                firstParam(componentNo),
                                MAX(lastParam(componentNo),normParam(componentNo)) 
                                 - firstParam(componentNo) + 1);          
                        break;
  	        case  7: /* cosmo */
                        outToTcl = realloc(outToTcl,(41+1)*sizeof(char));
			sprintf(outToTcl,"%13.5g %13.5g %13.5g",getH0(),getq0(),getl0());
			break;
                case  8: /* datagrp */
                        sprintf(outToTcl,"%8d", numberOfDataGroups());
                        break;
                case  9: /* datasets */
                        sprintf(outToTcl,"%8d", numberOfDataSets());
                        break;
                case 10: /* dof - degrees of freedom */
                        sprintf(outToTcl,"%8d",numberOfDataBins() - numberOfVariableParams());
                        break;
                /*case 11,12: see case 2 */
                /* options that require parameter number */
                case 13: /* error */
                case 27: /* param */
                case 29: /* pinfo */
                case 30: /* plink */
                        if (numberOfParams() > 1)
                        {
                                if ( objc < 3)
                                {
                                        char* errMessage = NULL;
                                        errMessage = (char *)malloc(64*sizeof(char));
                                        sprintf(errMessage,
                                         "Parameter number argument required for option %s",chosen);
                                        XWRITE(errMessage,10);
                                        free(errMessage);
                                        retStatus = TCL_ERROR;
                                        break;
                                         
                                }
                                next = Tcl_GetStringFromObj(objv[2],&arglen);
                                if (!isInteger(next))
                                {
                                        char errString[] = 
                                           "Argument for parameter number is not an integer";
                                        XWRITE(errString, 10);
                                        retStatus = TCL_ERROR;
                                        break;
                                }
                                sscanf(next,"%d",&paramNo);
                                if ( paramNo <= 0 || paramNo > numberOfParams())
                                {
                                        char* errMessage = NULL;
                                        errMessage = (char *)malloc(64*sizeof(char));
                                        sprintf(errMessage,
                                         "Parameter number p must satisfy 1 < p < %d",
                                         numberOfParams());
                                        XWRITE(errMessage,10);
                                        free(errMessage);
                                        retStatus = TCL_ERROR;
                                        break;
                                }
                        }
                        if (selected ==  13) /* error */
                        {
                                outToTcl = realloc(outToTcl,(48+1)*sizeof(char));
                                sprintf(outToTcl," %13.4g%13.4g %20s", paramValue(paramNo,"u"),
                                        paramValue(paramNo,"w"),confidenceIntervalError(paramNo));
                        }
                        if (selected == 27) /* param */
                        {
                                /* length of output is 7 - 14 character fields = 98 characters */
                                outToTcl = realloc(outToTcl,(98+1)*sizeof(char));
                                tmp = (char *)malloc((14+1)*sizeof(char));
                                sprintf(tmp," %13.4g",paramValue(paramNo,"v"));
                                strcpy(outToTcl,tmp);
                                sprintf(tmp," %13.4g",paramValue(paramNo,"d"));
                                strcat(outToTcl,tmp);
                                sprintf(tmp," %13.4g",paramValue(paramNo,"l"));
                                strcat(outToTcl,tmp);
                                sprintf(tmp," %13.4g",paramValue(paramNo,"b"));
                                strcat(outToTcl,tmp);
                                sprintf(tmp," %13.4g",paramValue(paramNo,"t"));
                                strcat(outToTcl,tmp);
                                sprintf(tmp," %13.4g",paramValue(paramNo,"h"));
                                strcat(outToTcl,tmp);
                                sprintf(tmp," %13.4g",paramValue(paramNo,"s"));
                                strcat(outToTcl,tmp);
                                free(tmp);
                        }
                        if (selected == 29) /* pinfo */
                        {
                                outToTcl = realloc(outToTcl,(18+1)*sizeof(char));
                                sprintf(outToTcl," %8s %8s",paramName(paramNo),paramUnit(paramNo));
                        }
                        if (selected == 30) /* plink */
                        {
                                if (paramIsLinked(paramNo))
                                {
                                        outToTcl = realloc(outToTcl,(34+1)*sizeof(char));
                                        sprintf(outToTcl," T %5d%13.4g%13.4g",
                                                linkParamNum(paramNo),addConToLink(paramNo),
                                                mulConToLink(paramNo));
                                }
                                else
                                {
                                        outToTcl = realloc(outToTcl,(1+1)*sizeof(char));
                                        sprintf(outToTcl,"F");
                                }
                                       
                        }
                        break;
                /*case 15,16: see case 2 */
                case 17: /* F statistic */
                        outToTcl = realloc(outToTcl,(29+1)*sizeof(char));
                        sprintf(outToTcl,"%14.6g %14.6g", getFStatisticValue(), getFStatisticProb());
                        break;
                case 18: /* genpop */
                        
                        next = Tcl_GetStringFromObj(objv[2],&arglen);
                        if (next)
                        {
                                if (isNumber(next))
                                {
                                        sscanf(next,"%d",&paramNo);
                                        if (paramNo <= 0 || paramNo > numberOfParams())
                                        {
                                                char* errMessage = NULL;
                                                errMessage = (char *)malloc(64*sizeof(char));
                                                sprintf(errMessage,
                                                  "Parameter number p must satisfy 1 < p < %d",
                                                  numberOfParams());
                                                XWRITE(errMessage,10);
                                                free(errMessage);
                                                retStatus = TCL_ERROR;
                                                break;
                                        }
                                }
                        }
                        else paramRequested = 0;

                        outToTcl = realloc(outToTcl,(14*genPopNum()+1)*sizeof(char));
                        tmp = (char *)malloc((14+1)*sizeof(char));
                        strcpy(outToTcl,"");
                        
                        if (paramRequested)
                        {
                                int fitParam = 0;
                                fitParam = fitParamNum(paramNo);
                                for (i = 0; i < genPopNum(); i++)
                                {
                                        sprintf(tmp," %13.5g",genFitPar(fitParam,i));
                                        strcat(outToTcl,tmp);
                                }
                                free(tmp);
                        }
                        else
                        {
                                for (i = 0; i < genPopNum(); i++)
                                {
                                        sprintf(tmp," %13.5g",genFitness(i));
                                        strcat(outToTcl,tmp);     
                                }        
                        }
                        break;
                case 19: /* goodness */
                        outToTcl = realloc(outToTcl,(14+1)*sizeof(char));
                        sprintf(outToTcl,"%14.6g", getGoodness());
                        break;
                case 20: /* idline */
                        if ( objc < 4)
                        {
                                char* errMessage = NULL;
                                errMessage = (char *)malloc(64*sizeof(char));
                                sprintf(errMessage,
                                 "Option %s requires 2 arguments: idline energy del(energy)",
                                   option[11]);
                                XWRITE(errMessage,10);
                                free(errMessage);
                                retStatus = TCL_ERROR;
                                break;
                                         
                        }
                        next = Tcl_GetStringFromObj(objv[2],&arglen);
                        
                        if ( !isNumber(next) )
                        {
                                char errString[] =  "Argument for energy is not a number";
                                XWRITE(errString, 10);
                                retStatus = TCL_ERROR;
                                break;
                        }
                        sscanf(next,"%f",&energy);
                        next = Tcl_GetStringFromObj(objv[3],&arglen);
                        if ( !isNumber(next) )
                        {
                                char errString[] =  "Argument for del(energy) is not a number";
                                XWRITE(errString, 10);
                                retStatus = TCL_ERROR;
                                break;
                        }
                        sscanf(next,"%f",&deltaEnergy);
                        tmp1 = (char *)malloc(256*sizeof(char));
                        tmp2 = (char *)malloc(256*sizeof(char));
                        getLineID(energy,deltaEnergy,1,tmp1,&status);
                        getLineID(energy,deltaEnergy,2,tmp2,&status);                        
                        outToTcl = realloc(outToTcl,(strlen(tmp1) + strlen(tmp2) + 1)*sizeof(char));
                        strcat(strcpy(outToTcl,tmp1),tmp2);
                        free(tmp1);
                        free(tmp2);
                        break;
                case 22: /* model  */
                        tmp = writeModel(25);
                        outToTcl = realloc(outToTcl,(strlen(tmp) + 1)*sizeof(char));
                        strcpy(outToTcl,tmp);
                        break;
                case 23: /* modcomp */
                        sprintf(outToTcl,"%8d",numberOfComponents());
                        break;
                case 24: /* modpar */
                        sprintf(outToTcl,"%8d",numberOfParams());
                        break;
                /* case  25, 26: see case 2 */
                /* case  27: see case 12 */
                /* case  28: see case 2 */
                /* cases 29,30: see case 12 */
                case 31: /* plot */
                        posn  = 5;
                        energies = udmArray("e");
                        components = udmArray("p");
                        flux     = udmArray("f");
                        group    = udmArray("g");
                        beginChan = udmArray("i");
                        endChan   = udmArray("l");
                        responseElts = udmArray("r");
                        effArea  = udmArray("a");
			tmp = (char *)malloc(250000*sizeof(char));
                        
                        plotit(inString, &posn, energies, components, flux, group, beginChan, 
                                endChan, responseElts, effArea, interactive, plot, tmp);
                        outToTcl = realloc(outToTcl, (strlen(tmp) + 1)*sizeof(char));
                        strcpy(outToTcl,tmp);
			free(tmp);
                        break;
                case 32: /* plotgrp */
                        sprintf(outToTcl,"%8d",numberOfPlotGroups());
                        break;
                case 33: /* query status */
                        outToTcl = realloc(outToTcl,(3+1)*sizeof(char));
		        if (askQuestion()) 
                        {
                           strcpy(outToTcl,"on");
			} else {
			   if (getAnswer())
			   {
			      strcpy(outToTcl,"yes");
			   } else {
			      strcpy(outToTcl,"no");
			   }
			}
                        break;
                /*case 34,35: see case 2 */
   	        case 36: /* simpars */
			nparams = numberOfParams();
  		        values = (double *)malloc(sizeof(double)*nparams);
			tmp = (char *)malloc((14+1)*sizeof(char));
                        strcpy(outToTcl,"");
                        outToTcl = realloc(outToTcl,sizeof(char)*(14*nparams+1));
			getSimulatedParamValues(nparams, values);
			for (i = 0; i < nparams; i++)
			  {
			    sprintf(tmp," %13.6g", *(values+i));
			    strcat(outToTcl,tmp);
			  }
			free(tmp);
			free(values);
			break;
                case 37: /* solab */
                        tmp = (char *)malloc(13+1*sizeof(char));
                        strcpy(outToTcl,"");
                        for (i = 0; i < 18; ++i)
                        {
                                outToTcl = realloc(outToTcl,((i+1)*13 + 1)*sizeof(char));
                                sprintf(tmp," %12.4g",abundance(celt[i]));
                                strcat(outToTcl,tmp);
                        }
                        free(tmp);
                        break;
                case 38: /* stat - ftstat is input parameter. */
                        outToTcl = realloc(outToTcl,(14+1)*sizeof(char));
                        sprintf(outToTcl,"%14.6g",getFitStatistic());
                        break;
                case 39: /* varpar */
                        sprintf(outToTcl,"%8d",numberOfVariableParams());   
                        break;     
        }
        free(inString);
/*        free(outToTcl); */
        free(chosen);
        if (retStatus == TCL_OK) Tcl_SetVar(interp,"xspec_tclout",outToTcl,0);  
        return retStatus;



}