#include <UserInterface/xstcl.h>
#include <math.h>

#include <XSUser/Global/Global.h>
#include <XSUser/Global/XSGlobal.h>
#include <XSFit/Fit/Fit.h>
#include <XSFit/Fit/FitMethod.h>
#include <XSFit/Fit/StatMethod.h>
#include <XSFit/MCMC/AsciiChain.h>
#include <XSFit/MCMC/FITSChain.h>
#include <XSFit/MCMC/ChainManager.h>
#include <XSFit/MCMC/MarginGrid.h>
#include <XSFit/Randomizer/RandomizerBase.h>
#include <XSFunctions/FunctionUtility.h>
#include <XSModel/GlobalContainer/DataContainer.h>
#include <XSModel/GlobalContainer/ModelContainer.h>
#include <XSModel/GlobalContainer/ModelTypes.h>
#include <XSModel/GlobalContainer/ResponseContainer.h>
#include <XSModel/Data/SpectralData.h>
#include <XSModel/Data/BackCorr/Background.h>
#include <XSModel/Data/Detector/Response.h>
#include <XSModel/Data/Detector/RealResponse.h>
#include <XSModel/Data/Detector/MultiResponse.h>
#include <XSModel/Model/Model.h>
#include <XSModel/Model/EmissionLines/LineList.h>
#include <XSModel/Model/Component/Component.h>
#include <XSModel/Parameter/ModParam.h>
#include <XSPlot/Plot/Plot.h>
#include <XSUtil/Numerics/Numerics.h>
#include <XSUtil/Utils/XSutility.h>
#include <XSUtil/Utils/XSstream.h>
#include <XSUtil/Parse/XSparse.h>

#include <XSstreams.h>
#include <XSContainer.h>
#include <XSUser/Handler/HandlerUtils.h>
#include <memory>
#include <sstream>
  
int
XSGlobal::xsTclout(ClientData cdata,Tcl_Interp* interp,int objc, Tcl_Obj* CONST objv[] )
{
   using namespace std;
   using namespace XSContainer;        

   enum SUBCOMS {HELP, AREASCAL, ARF, BACKGRND, BACKSCAL, CHAIN, 
        CHATTER, COMPINFO, COSMO, COVARIANCE, DATAGRP, DATASETS, DOF, ENERGIES, 
        EQWIDTH, ERROR, EXPOS, FILENAME, FLUX, FTEST, GENPOP, GOODNESS, 
        IDLINE, LUMIN, MARGIN, MODEL, MODCOMP, MODPAR, MODVAL, NOTICED, PARAM, 
        PEAKRSID, PINFO, PLINK, PLOT, PLOTGRP, QUERY, RATE, RESPONSE, 
        SIMPARS, SOLAB, STAT, STATMETHOD, STEPPAR, VARPAR, WEIGHT, XFLT};
   static map<string, size_t> option;
   const string cmd (static_cast<char *>(cdata));
   
   if (!option.size())
   {
      option["?"]          = HELP;
      option["areascal"]   = AREASCAL;
      option["arf"]        = ARF;
      option["backgrnd"]   = BACKGRND;
      option["backscal"]   = BACKSCAL;
      option["chain"]      = CHAIN;
      option["chatter"]    = CHATTER;
      option["compinfo"]   = COMPINFO;
      option["cosmo"]      = COSMO;
      option["covariance"] = COVARIANCE;
      option["datagrp"]    = DATAGRP;
      option["datasets"]   = DATASETS;
      option["dof"]        = DOF;
      option["energies"]   = ENERGIES;
      option["eqwidth"]    = EQWIDTH;
      option["error"]      = ERROR;
      option["expos"]      = EXPOS;
      option["filename"]   = FILENAME;
      option["flux"]       = FLUX;
      option["ftest"]      = FTEST;
      option["genpop"]     = GENPOP;
      option["goodness"]   = GOODNESS;
      option["idline"]     = IDLINE;
      option["lumin"]      = LUMIN;
      option["margin"]     = MARGIN;     
      option["model"]      = MODEL;
      option["modcomp"]    = MODCOMP;
      option["modpar"]     = MODPAR;
      option["modval"]     = MODVAL;
      option["noticed"]    = NOTICED;
      option["param"]      = PARAM;
      option["peakrsid"]   = PEAKRSID;
      option["pinfo"]      = PINFO;
      option["plink"]      = PLINK;
      option["plot"]       = PLOT;
      option["plotgrp"]    = PLOTGRP;
      option["query"]      = QUERY;
      option["rate"]       = RATE;
      option["response"]   = RESPONSE;
      option["simpars"]    = SIMPARS;
      option["solab"]      = SOLAB;
      option["stat"]       = STAT;
      option["statmethod"] = STATMETHOD;
      option["steppar"]    = STEPPAR;
      option["varpar"]     = VARPAR;
      option["weight"]     = WEIGHT;
      option["xflt"]       = XFLT; 
   }
   
   if (objc == 1)  
   {
      XSutility::printValidOptions(tcout, cmd, option);
      return TCL_OK;
   }
   
   string inputOpt(Tcl_GetString(objv[1]));
   map<string,size_t>::const_iterator itOpt(option.lower_bound(inputOpt));
   map<string,size_t>::const_iterator itOptEnd(option.end());
   
   if (itOpt != itOptEnd && itOpt->first.find(inputOpt) == 0)
   {
      const size_t selected = itOpt->second;
      const string& opt = itOpt->first;
      size_t specNum=0;
      string outToTcl;
      std::ostringstream ss;
      char XSPEC_TCLOUT[] = "xspec_tclout";

      switch (selected)
      {
         case HELP:  // ?
            XSutility::printValidOptions(tcout, cmd, option);
            break;
         case AREASCAL: // areascal n s|b (default to s)
         case BACKSCAL: // backscal n s|b
         case EXPOS:    // exposure n s|b
            if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
            {
               return TCL_ERROR;
            }
            else
            {
               const SpectralData* sd = datasets->lookup(specNum);
               bool isBackground = false;
               const RealArray* pScaleArray = 0;
               Real exposTime = -999.0;
               if (objc > 3)
               {
                  string sOrB(Tcl_GetString(objv[3]));
                  if (sOrB[0] == 'b' || sOrB[0] == 'B')
                        isBackground = true;
               }
               if (isBackground)
               {
                  const Background* back = sd->background();
                  if (!back)
                  {
                     tcerr << "No background data exists for spectrum " 
                           << specNum << std::endl;
                     return TCL_ERROR;
                  }
                  const SpectralData* bsd = back->data();
                  if (selected == AREASCAL|| selected == BACKSCAL)
                  {
                     pScaleArray = (selected == AREASCAL) ?
                        &bsd->areaScale() : &bsd->backgroundScale();
                  }
                  else
                  {
                     exposTime = bsd->exposureTime();
                  }
               }
               else
               {
                  if (selected == AREASCAL|| selected == BACKSCAL)
                  {
                     pScaleArray = (selected == AREASCAL) ?
                        &sd->areaScale() : &sd->backgroundScale();
                  }
                  else
                  {
                     exposTime = sd->exposureTime();
                  }
               }
               
               if (exposTime >= .0)
               {
                  ss << exposTime;
               }
               else
               {
                  size_t sz = pScaleArray->size();
                  for (size_t i=0; i<sz; ++i)
                  {
                     ss << (*pScaleArray)[i] << "   ";
                  }
               }
            }            
            break;   
         case ARF: // arf n
            if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
            {
               return TCL_ERROR;               
            }
            else
            {
               SpectralData* sd = datasets->lookup(specNum);
               string arfNames("");
               const std::vector<Response*>& resps = sd->detector();
               for (size_t i=0; i<resps.size(); ++i)
               {
                  const RealResponse* rrsp = dynamic_cast<RealResponse*>(resps[i]);
                  if (rrsp)
                  {
                     arfNames += rrsp->arfName();
                     arfNames += " ";
                  }
                  else
                  {
                     const MultiResponse* mrsp = dynamic_cast<MultiResponse*>(resps[i]);
                     if (mrsp)
                     {
                        const StringArray& arfs = mrsp->arfNames();
                        for (size_t j=0; j<arfs.size(); ++j)
                        {
                           arfNames += arfs[j];
                           arfNames += " ";
                        }
                     }
                  }
               }
               ss << arfNames;
            }
            break;
	 case BACKGRND: // backgrnd n
            if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
            {
               return TCL_ERROR;
            }
            else
	    {
	       SpectralData* sd = datasets->lookup(specNum);
	       const string& backFile = sd->backgroundFile();
               ss << backFile;
	    }
	    break;
         case CHAIN: // last|proposal|stat
            {
               const string chainLast("last");
               const string chainProposal("proposal");
               const string chainStat("stat");
               if (objc > 2)
               {
                  string secondOpt = XSutility::lowerCase(string(Tcl_GetString(objv[2])));
                  if (chainLast.find(secondOpt)==0)
                  {
                     if (objc > 3)
                     {
                        const string fileName = string(Tcl_GetString(objv[3]));
                        // First see if this file happens to be a loaded chain.
                        // If not, attempt to create a temporary new chain
                        // from it.
                        const Chain* chain = fit->chainManager()->getChain(fileName);
                        std::auto_ptr<Chain> apChain(0);
                        if (!chain)
                        {
                           // Create a new Chain.  We don't know the file format 
                           // or even if it exists.
                           std::auto_ptr<ChainIO> apIO(0);
                           try
                           {
                              apIO.reset(new AsciiChain(fileName));
                           }
                           catch (ChainIO::ChainIOError&)
                           {
                              // If in here, file cannot even be opened.
                              return TCL_ERROR;
                           }
                           catch (YellowAlert&)
                           {
                              // File exists, but not in the proper text format.
                              try
                              {
                                 apIO.reset(new FITSChain(fileName));
                              }
                              catch (YellowAlert&)
                              {
                                 tcerr << "***Error: Chain file has an unrecognized format."
                                     << std::endl;
                                 return TCL_ERROR;
                              }
                           }
                           // File exists and we now know which format to try.
                           try
                           {
                              apChain.reset(new Chain(apIO.get(), fileName));
                              apIO.release();
                           }
                           catch (YellowAlert&)
                           {
                              return TCL_ERROR;
                           }
                           // apChain will still own the memory.
                           chain = apChain.get();
                        }
                        try
                        {
                           // A chain object can't exist with width = 0, so
                           // this is safe.
                           RealArray lastPoint(.0, chain->width()-1);
                           chain->openForReadPoints();
                           chain->readPointFromLine(chain->length()-1, lastPoint);
                           chain->closeFile();
                           for (size_t i=0; i<lastPoint.size(); ++i)
                              ss << lastPoint[i] <<"  ";
                        }
                        catch (YellowAlert&)
                        {
                           if (chain)
                           {
                              // If this is a temporary chain, its file would
                              // automatically close as soon as auto_ptr goes
                              // out of scope.  But if chainManager owns it,
                              // file will be left open if something threw after
                              // openForReadPoints call.
                              // This call does nothing if file is already closed.
                              chain->closeFile();
                           }
                           return TCL_ERROR;
                        }
                     }
                     else
                     {
                        tcerr << "Tclout chain last also requires a file name argument."
                           << std::endl;
                        return TCL_ERROR;
                     }
                  }
                  else if (chainProposal.find(secondOpt)==0) // distribution|matrix
                  {
                     const string chainDist("distribution");
                     const string chainMatrix("matrix");
                     if (objc > 3)
                     {
                        string infoSpec = XSutility::lowerCase(string(Tcl_GetString(objv[3])));
                        if (chainDist.find(infoSpec)==0)
                        {
                           ss << fit->chainManager()->chainProposal()->name();
                        }
                        else if (chainMatrix.find(infoSpec)==0)
                        {
                           const RealArray* pCovarMatrix = 
                              fit->chainManager()->chainProposal()->covariance(fit);
                           if (pCovarMatrix)
                           {
                              const RealArray& covarMatrix = *pCovarMatrix;
                              // For consistency with chain prop covariance input,
                              // only print out the lower half and diagonal.
                              // We're going to ASSUME this is a square matrix.
                              const size_t nPar = static_cast<size_t>
                                      (sqrt(static_cast<Real>(covarMatrix.size())));
                              for (size_t i=0; i<nPar; ++i)
                              {
                                 for (size_t j=0; j<=i; ++j)
                                 {
                                    ss << covarMatrix[i*nPar+j] << " ";
                                 }
                              }
                           }
                        }
                        else
                        {
                           tcerr << "Valid chain proposal args: distribution|matrix"
                              << std::endl;
                           return TCL_ERROR;
                        } 
                     }
                     else
                     {
                        tcerr << "Valid chain proposal args: distribution|matrix"
                           <<std::endl;
                        return TCL_ERROR;
                     } 
                  } // end if chain proposal 
                  else if (chainStat.find(secondOpt)==0)
                  {
                     ss << fit->chainManager()->lastStatCalc();
                  } // end if chain stat
                  else
                  {
                     tcerr << "Valid chain options: last | proposal | stat" << std::endl;
                     return TCL_ERROR;
                  }
               } // end if objc > 2
               else
               {
                  tcerr << "Valid chain options: last | proposal | stat" << std::endl;
                  return TCL_ERROR;
               }
            }          
            break;
	 case CHATTER: // chatter
	    {
               ss << tpout.consoleChatterLevel() << " " 
                  << tpout.logChatterLevel();
	    }
	    break;
         case COMPINFO: // compinfo <modelName>:n <dataGroup n>
            {
               string modName;
               size_t idx = string::npos;
               const Component* comp = 0;
               if (HandlerUtils::verifyStringInt<Component>(objc, objv, 3, modName, idx))
               {
                  comp = HandlerUtils::lookupStringIntObj<Component>(modName, idx);
               }
               if (!comp)
               {
                  return TCL_ERROR;
               }
               else
               {
                  size_t dgNum = 1;
                  size_t groupPos = 0;
                  if (objc > 3)
                  {
                     size_t tmpNum = XSutility::isInteger(string(Tcl_GetString(objv[3])));
                     if (tmpNum != string::npos)
                     {
                        dgNum = tmpNum;
                     }                     
                  }
                  const Model* mod = HandlerUtils::getModFromName(modName);
                  // If the earlier call to retrieve the component
                  // succeeded, there should be no way for this to
                  // fail.  Still, just to be safe.....
                  if (!mod)
                  {
                     return TCL_ERROR;
                  }
                  if (!mod->isActive())
                  {
                     if (dgNum > 1)
                     {
                        tcerr << "Invalid data group number specified."
                           <<"\n Model is currently not active" << std::endl;
                        return TCL_ERROR;
                     }
                     groupPos = 1;
                  }
                  else
                  {
                     std::map<size_t,std::map<size_t,size_t> >::const_iterator 
                           itSource = datasets->sourceToDgs().find(mod->sourceNumber());
                     if (itSource != datasets->sourceToDgs().end())
                     {
                        std::map<size_t,size_t>::const_iterator itGroup = 
                                itSource->second.find(dgNum);
                        if (itGroup != itSource->second.end())
                        {
                           groupPos = itGroup->second;
                        }
                     }
                     if (!groupPos)
                     {
                        tcerr << "Invalid data group number specified for model." 
                           << std::endl;
                        return TCL_ERROR;
                     }
                  }
                  size_t nTotPar = mod->numberOfParameters();
                  const std::vector<Parameter*>& compPars = comp->parameterSet();
                  size_t nPar = compPars.size();
                  size_t firstPar = 0;
                  if (nPar)
                  {
                     firstPar = compPars[0]->index();
                     firstPar += (groupPos-1)*nTotPar;
                  }
                  ss << comp->name() << " " << firstPar << " " <<nPar;
               }
            }
            break;
         case COSMO:  // cosmo
            ss << models->cosmo().H0 << "   " << models->cosmo().q0
               << "   " << models->cosmo().lambda0;
            break;
         case COVARIANCE:  // covariance [m, n]
            {
               const RealArray& covar = fit->fitMethod()->covariance();
               size_t sz = covar.size();
               if (sz)
               {
                  // Just assumes covar is NxN or empty.
                  size_t N = static_cast<size_t>(sqrt(static_cast<Real>(sz)));
                  if (objc > 2)
                  {
                     StringArray inArgs, outArgs;
                     IntegerArray dummy;
                     for (int i=2; i<objc; ++i)
                     {
                        inArgs.push_back(Tcl_GetString(objv[i]));
                     }
                     XSparse::collectParams(inArgs, dummy, outArgs);
                     if (outArgs.size() < 2)
                     {
                       tcerr << "2 indices are required to specify covariance element (m, n)"<<endl;
                       return TCL_ERROR;
                     }
                     size_t m = XSutility::isInteger(outArgs[0]);
                     size_t n = XSutility::isInteger(outArgs[1]);
                     if (m == string::npos || n == string::npos)
                     {
                        tcerr << "Improper covariance index specifier." << endl;
                        return TCL_ERROR;
                     }
                     if (m == 0 || n == 0 || m > N || n > N)
                     {
                        tcerr << "Covariance index is out of range."
                          << "  Valid range is 1 - " << N << endl;
                        return TCL_ERROR;
                     }
                     ss << covar[(m-1)*N + n-1];
                  }
                  else
                  {
                     for (size_t i=0; i<N; ++i)
                     {
                        for (size_t j=0; j<=i; ++j)
                        {
                           ss << covar[i*N+j] << " ";
                        }
                     }
                  }
               }
               else
               {
                  tcerr << "No covariance is currently calculated." << endl;
                  return TCL_ERROR;
               }
            }
            break;
	 case DATAGRP: // datagrp [n]
            if (objc < 3)
            {
               ss << datasets->numberOfGroups();
            }
            else
            {
               if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
               {
                  return TCL_ERROR;
               }
               else
               {
                  size_t dummyRow=0;
                  const DataSet* ds = datasets->dataSetLookup(specNum, dummyRow);
                  ss << ds->dataGroup();
               }
            }
            break;
	 case DATASETS: // datasets	    
            ss << datasets->numberOfSpectra();
            break;
         case DOF: // dof
            {
               int degrees(0);
               size_t bins(0);
               fit->statMethod()->degreesOfFreedom(fit,degrees,bins);
               ss << degrees << " " << bins;
            }
            break;
         case ENERGIES: // energies n
            if (models->autonomousEnergy().size())
            {
               const RealArray& engs = models->autonomousEnergy();
               for (size_t i=0; i<engs.size(); ++i)
               {
                  ss << engs[i] << "  ";
               }
            }
            else
            {
               if ((!datasets->numberOfSpectra() && objc == 2) || (objc > 2
                           && string(Tcl_GetString(objv[2]))[0] == '0' ))
               {
                  // use singleton dummy response
                  const Response* dummyrsp = responses->responseList(DUMMY_RSP,1);
                  if (dummyrsp)
                  {
                     const RealArray& engs = dummyrsp->energies();
                     for (size_t i=0; i<engs.size(); ++i)
                     {
                        ss << engs[i] << "   ";
                     }
                  }               
               }
               else
               {
                  if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
                  {
                     return TCL_ERROR;
                  }
                  else
                  {
                     // Hardcoding this to 1 until multiple sources are 
                     // implemented for this option.
                     const size_t sourceNum = 1;
                     // If extended energies are in use, energy arrays must
                     // come from the model (if it exists) since responses 
                     // don't know about extensions.  
                     const SpectralData* sd = datasets->lookup(specNum);
                     const Response* rsp = sd->detector(sourceNum-1);
                     if (rsp)
                     {
                        const RealArray *engs=0;
                        RealArray ext;
                        if (models->isExtended())
                        {
                           const Model* mod = models->lookup(rsp);
                           if (mod)
                           {
                               ArrayContainer::const_iterator itEngs = 
                                        mod->energy().find(specNum);
                               engs = &(itEngs->second);
                           }
                           else
                           {
                              ext.resize(rsp->energies().size());
                              ext = rsp->energies();
                              Model::makeExtendArray(models->extendedEnergy(), ext);
                              engs = &ext;                                        
                           }
                        }
                        else
                        {
                           engs = &rsp->energies();
                        }
                        for (size_t i=0; i<engs->size(); ++i)
                        {
                           ss << (*engs)[i] << "   ";
                        }
                     }
                     else
                     {
                        tcerr << "No response associated with spectrum "
                           << specNum << " source 1." << std::endl;
                        return TCL_ERROR;
                     }
                  }
               }
            }            
            break;   
         case EQWIDTH: // eqwidth n
            if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
            {
               return TCL_ERROR;               
            }
            else
            {
               const SpectralData* sd = datasets->lookup(specNum);
               const SpectralData::FluxCalc& eqWidth = sd->lastEqWidthCalc();
               ss << eqWidth.value << " " << eqWidth.errLow << " " 
                  << eqWidth.errHigh;
            }
            break;
         case ERROR: // error <modelName:>n errCode
            {
               const Parameter* param = HandlerUtils::getFromStringInt<Parameter>(objc, objv);
               if (!param)
               {
                  return TCL_ERROR;
               }
               try
               {
                  const ModParam* mPar = dynamic_cast<const ModParam*>(param);
                  ss << mPar->emn() << " " << mPar->epe() << "  "
                        << mPar->lastErrorStatus();
               }
               catch (std::bad_cast)
               {
                  tcerr << " Parameter is not a fit parameter." << std::endl;
                  return TCL_ERROR;
               }
               
            }
            break;
         case FILENAME: // filename n
            if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
            {
               return TCL_ERROR;
            }
            else
            {
               size_t dum=0;
               const DataSet* ds = datasets->dataSetLookup(specNum, dum);
               string fileName = ds->getFullPathName();
               ss << fileName;
            }
            break;
         case FLUX: // flux [n]
	     if(datasets->numberOfSpectra())
	     {
                if (datasets->numberOfSpectra() > 1 && objc == 2)
                {
                   tcerr << "***XSPEC Error: More than 1 spectrum is currently loaded."
                     <<"\n   Must specify spectrum number n: \"tclout flux <n>\"" <<std::endl;
                   return TCL_ERROR;
                }  
		if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
                {
                   return TCL_ERROR;
                }
		const SpectralData* sd = datasets->lookup(specNum);
		size_t nFlux = sd->lastModelFluxCalc().size();

		for (size_t i=0; i<nFlux; ++i)
		{
		    const SpectralData::FluxCalc& flux = sd->lastModelFluxCalc(i);
		    ss << flux.value <<" "<< flux.errLow <<" "<< flux.errHigh <<"  "
                       << flux.photonValue <<" "<< flux.photonLow <<" "
                       << flux.photonHigh <<" ";
		}
		if (!nFlux)
		    ss << "0.  0.  0.  0.  0.  0.";
	     }
	     else
	     {
		 //get vector of Model pointers
		 std::vector<Model*> mods = XSutility::values<
		     std::map<string, Model*>, Model*>(models->modelSet().begin(),
						       models->modelSet().end());

		 if(mods.size())
		 {
		     //sort based on sourcenumber using binary predicate CompareSource()
		     std::sort(mods.begin(), mods.end(), HandlerUtils::CompareSource());

		     std::vector<Model*>::iterator
			 begMods = mods.begin(), endMods = mods.end();

		     while(begMods != endMods)
		     {
			 if(models->activeModelNames((*begMods)->name()))
			 {
			     Model* m = *begMods;
			 
			     const SpectralData::FluxCalc& lumin = m->lastModelFluxCalc();
			     ss << lumin.value <<" "<< lumin.errLow <<" "<< lumin.errHigh<<"  "
                               << lumin.photonValue <<" "<< lumin.photonLow <<" "
                               << lumin.photonHigh <<" ";
			 }
			 ++begMods;
		     }
		 }
		 else
		 {
		     tcerr << "No spectra are loaded: Unable to process  " << opt 
			   << " option" << std::endl;
		     
		     return TCL_ERROR;
		 }
	     }
	     break;
         case FTEST: // ftest
            ss << Fit::lastFtest();
            break;
         case GOODNESS:
            ss << Fit::lastGoodness();
            break;
         case IDLINE: // idline e d
            if (objc < 4)
            {
               tcerr << "Option requires 2 arguments: idline energy del(energy)"
                     << std::endl;
               return TCL_ERROR; 
            }
            else
            {
               RealArray input(.0,2);
               for (size_t i=0; i<2; ++i)
               {
                  string arg(Tcl_GetString(objv[i+2]));
                  if (!XSutility::isReal(arg, input[i]))
                  {
                     tcerr << "Argument is not a number" << std::endl;
                     return TCL_ERROR;
                  }
               }
               Real lowE = input[0] - input[1];
               Real highE = input[0] + input[1];
               LineList *lineList = 0;
               try
               {
                  lineList = LineList::get("bearden");
                  if (lineList)
                  {
                     std::ostringstream ssTmp;
                     lineList->initialize(lowE, highE, .0, .0, false);
                     // This can throw:
                     lineList->showList(ssTmp);
                     outToTcl = ssTmp.str();                     
                  }
                  else
                  {
                     tcerr << "***Warning: \"bearden\" is an unrecognized list name."
                            << std::endl;
                  }
                  lineList = LineList::get("mekal");
                  if (lineList)
                  {
                     std::ostringstream ssTmp;
                     lineList->initialize(lowE, highE, .0, .0, false);
                     // This can throw:
                     lineList->showList(ssTmp);
                     outToTcl += ssTmp.str();                     
                  }
                  else
                  {
                     tcerr << "***Warning: \"mekal\" is an unrecognized list name."
                            << std::endl;
                  }
                  ss << outToTcl;
               }
               catch (YellowAlert&)
               {
                  delete lineList;
                  return TCL_ERROR;
               }
               delete lineList;
            }
            break;
         case LUMIN: // lumin n
	     if(datasets->numberOfSpectra())
	     {
		 HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum);
		 
		 const SpectralData* sd = datasets->lookup(specNum);
		 
		 size_t nFlux = sd->lastModelLuminCalc().size();
		 for (size_t i=0; i<nFlux; ++i)
		 {
		     const SpectralData::FluxCalc& lumin = sd->lastModelLuminCalc(i);
		     ss << lumin.value <<" "<< lumin.errLow <<" "<< lumin.errHigh<<"  "
                        << lumin.photonValue <<" "<< lumin.photonLow <<" "
                        << lumin.photonHigh <<" ";
		 }
		 if (!nFlux)
		 {
		     ss << "0.  0.  0.  0.  0.  0.";
		 }
	     }
	     else
	     {
		 std::vector<Model*> mods = XSutility::values<
		     std::map<string, Model*>, Model*>(models->modelSet().begin(), 
						       models->modelSet().end());

		 if(mods.size())
		 {
		     //sort based on sourcenumber using binary predicate CompareSource()
		     std::sort(mods.begin(), mods.end(), HandlerUtils::CompareSource());

		     std::vector<Model*>::iterator
			 begMods = mods.begin(), endMods = mods.end();

		     while(begMods != endMods)
		     {
			 if(models->activeModelNames((*begMods)->name()))
			 {
			     Model* m = *begMods;
			 
			     const SpectralData::FluxCalc& lumin = m->lastModelLuminCalc();
			     ss << lumin.value <<" "<< lumin.errLow <<" "<< lumin.errHigh<<"  "
                                << lumin.photonValue <<" "<< lumin.photonLow 
                                <<" " << lumin.photonHigh << " ";
			 }
			 ++begMods;
		     }
		 }
		 else
		 {
		     tcerr << "No spectra are loaded: Unable to process  " << opt 
			   << " option" << std::endl;
		     
		     return TCL_ERROR;
		 }
	     }
	     break;
         case MODEL: //model
            if (!XSContainer::models->modelSet().empty())
            {
               const ModelMap& mods = XSContainer::models->modelSet();
               ModelMapConstIter mod = mods.begin();
               ModelMapConstIter modEnd = mods.end();
               while (mod != modEnd)
               {
                  if (mod->second->name() != Model::DEFAULT())
                  {
                     ss << mod->second->name() << ": ";
                  }
                  ss << mod->second->fullExpression() << std::endl;
                  ++mod;
               }
            }
            break;
         case MODCOMP: // modcomp name
            {
               const Model* mod = HandlerUtils::getModFromArg(objc, objv);
               if (!mod)
                  ss << 0;
               else
                  ss << mod->numberOfComponents();
            }
            break;
         case MODPAR:  // modpar name
            {
               const Model* mod = HandlerUtils::getModFromArg(objc, objv);
               if (!mod)
                  ss << 0;
               else
               {
                  size_t numPar = mod->numberOfParameters();
                  if (mod->isActive())
                  {
                     numPar *= datasets->getNumberOfGroupsForSource(mod->sourceNumber());
                  }
                  ss << numPar;
               }
            }
            break;
         case MODVAL:  // modval [<specNum> [<modName>]]
            if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
            {
               return TCL_ERROR;               
            }
            else
            {
               string modName(Model::DEFAULT());
               string errMsg("Default (unnamed) model");
               if (objc >= 4)
               {
                  modName = string(Tcl_GetString(objv[3]));
                  errMsg = string("Model ") + modName;
               }
               size_t dummyRow=0;
               size_t dgNum = datasets->dataSetLookup(specNum, dummyRow)->dataGroup();
               const Model* mod = models->lookup(modName, dgNum);
               if (!mod)
               {
                  tcerr << errMsg << " not found associated with spectrum " 
                        << specNum << endl;
                  return TCL_ERROR;
               }
               try 
               {
                  const RealArray& modFlux = mod->modelFlux(specNum);
                  size_t nE = modFlux.size();
                  for (size_t i=0; i<nE; ++i)
                  {
                     ss << modFlux[i] << "  ";
                  }
               }
               catch (YellowAlert&)
               {
                  tcerr << errMsg << " not found associated with spectrum " 
                        << specNum << endl;
                  return TCL_ERROR;
               }
            }
            break;
         case NOTICED: // noticed [<specNum>] | energy [<specNum>]
         {
            const string enString("energy");
            string enTest;
            if (objc > 2)
               enTest = XSutility::lowerCase(string(Tcl_GetString(objv[2])));
            if (objc > 2 && enString.find(enTest) == 0)
            {
               // output noticed energy ranges
               if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 4, specNum))
               {
                  return TCL_ERROR;               
               }
               else
               {
                  const SpectralData* sd = datasets->lookup(specNum);
                  const BoolArray& noticedChans = sd->noticedChannels();
                  std::vector<Response*>::const_iterator itDet = 
                        sd->detector().begin();
                  std::vector<Response*>::const_iterator itDetEnd = 
                        sd->detector().end();
                  const Response* firstResp=0;
                  while (!firstResp && itDet != itDetEnd)
                  {
                     if (*itDet)
                        firstResp = *itDet;
                     ++itDet;
                  }
                  if (!firstResp)
                  {
                     tcerr << "\nNo response associated with spectrum " << specNum
                        <<".  Cannot convert channels to energies." << std::endl;
                     return TCL_ERROR;               
                  }
                  const RealArray& eMin = firstResp->eboundsMin();
                  const RealArray& eMax = firstResp->eboundsMax();
                  const size_t nChans = sd->channels();
                  const size_t offset = sd->startChan()-sd->firstChan();
                  bool isPrevNoticed = false;
                  bool isWave = SpectralData::wavelengthSpace();
                  Real startRange=0.0, endRange=0.0;
                  for (size_t i=0; i<nChans; ++i)
                  {
                     if (noticedChans[i+offset])
                     {
                        if (!isPrevNoticed)
                        {
                           startRange = eMin[i];
                           if (isWave)
                           {
                              startRange = std::max(SMALL,startRange);
                              startRange = Numerics::KEVTOA/startRange;
                           }
                           ss << startRange <<"-";
                        }
                        endRange = eMax[i];
                        isPrevNoticed = true;
                     }
                     else
                     {
                        if (isPrevNoticed)
                        {
                           if (isWave)
                           {
                              endRange = std::max(SMALL,endRange);
                              endRange = Numerics::KEVTOA/endRange;
                           }
                           ss << endRange <<" ";
                        }
                        isPrevNoticed = false;
                     }
                  }
                  if (isPrevNoticed)
                  {
                     if (isWave)
                     {
                        endRange = std::max(SMALL,endRange);
                        endRange = Numerics::KEVTOA/endRange;
                     }
                     ss << endRange;
                  }
               }  
            }
            else
            {
               // output noticed channel numbers
               if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
               {
                  return TCL_ERROR;               
               }
               else
               {
                  std::valarray<size_t> indirectNotice;
                  const SpectralData* sd = datasets->lookup(specNum);
                  sd->buildIndirectNotice(indirectNotice);
                  size_t n = indirectNotice.size();
                  // Using a for loop rather than just valarray += 1 due to some
                  // unexplained compiler issue on Solaris.
                  for (size_t i=0; i<n; ++i)
                  {
                     indirectNotice[i] += 1;
                  }
                  ss << XSparse::ArrayToRange(indirectNotice);
               }
            }
         }
            break;
         case PARAM: // param modelName:n
            {
               const Parameter* param = HandlerUtils::getFromStringInt<Parameter>(objc, objv);
               if (!param)
               {
                  return TCL_ERROR;
               }
               const string delim(" ");
               const string fullString = param->stringVal();
               string valueString;
               // 6 values for ModParams, 1 for Switch and Scale.
               size_t nVals = 6;
               if (!dynamic_cast<const ModParam*>(param)) nVals = 1;
               string::size_type iPos = 0;
               for (size_t i=0; i<nVals; ++i)
               {
                  iPos = fullString.find_first_not_of(delim, iPos);
                  iPos = fullString.find_first_of(delim, iPos);
               }
               valueString = fullString.substr(0, iPos); 
               ss << valueString;
            }
            break;
         case PEAKRSID: // peakrsid n {lo,hi}
            if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
            {
               return TCL_ERROR;
            }
            else
            {
               try   // fit::peakResidual can throw 
               {
                  using namespace XSutility;
                  std::pair<Real,Real> range(0,0);       
                  if ( objc >= 5)
                  {
                     string argn(Tcl_GetString(objv[3]));
                     string argn2(Tcl_GetString(objv[4]));
                     if (!isReal(argn,range.first) || !isReal(argn2,range.second))
                     {
                          tcerr << " syntax error: arguments following spectrum number "
                                << " must be real energy range. \n";
                          return TCL_ERROR;
                     }
                  }               
                  std::pair<Real,Real> max(0,0);
                  std::pair<Real,Real> min(0,0);
                  fit->peakResidual(max,min,specNum,range);
                  ss << max.second <<"  "<<max.first<<"  "<<min.second<<"  "
                       <<min.first;
               }
               catch (YellowAlert&)
               {
                  return TCL_ERROR;
               }
            }
            
            break;
         case PINFO: // pinfo modelName:n
            {
               const Parameter* param = HandlerUtils::getFromStringInt<Parameter>(objc, objv);
               if (!param)
               {
                  return TCL_ERROR;
               }
               string parVal = param->stringVal();
               int c(0);
               int ignore(7);
               if ( !dynamic_cast<const ModParam*>(param)) ignore = 1;

               // ignore the first "ignore" arguments, 
               // 7 for ModParams and 1 for Switch and Scale types.
               while ( c++ < ignore) XSparse::returnDelimitedArgument(parVal," ");
               ss << parVal;
            }
            break;
         case PLINK: // plink modelName:n
            {
               const Parameter* par = HandlerUtils::getFromStringInt<Parameter>(objc, objv);
               if (!par)
               {
                  return TCL_ERROR;
               }
               if (par->isLinked() )
               {
                  const string link(par->parameterSetting());
                  ss << "T " << link;     
               }
               else
               {
                  ss << "F";
               }
            }
            break;
         case PLOT: // plot <option> <array> [<plot group n>]
         {
            map<string, Plot::SaveArrayOption> arrayOpt;
            arrayOpt["x"] = Plot::X_AR;
            arrayOpt["xerr"] = Plot::XERR_AR;
            arrayOpt["y"] = Plot::Y_AR;
            arrayOpt["yerr"] = Plot::YERR_AR;
            arrayOpt["model"] = Plot::MODEL_AR;
            static string savPlotCom = "data";
            static Plot::SaveArrayOption savArrayEnum = Plot::X_AR;
            static int savPlotGroup = 1;
            
            StringArray tclInArgs;
            StringArray params;
            IntegerArray iParams;
            for (int i=2; i<objc; ++i)
            {
               tclInArgs.push_back(Tcl_GetString(objv[i]));
            }
            XSparse::collectParams(tclInArgs, iParams, params);
            map<string, Plot::SaveArrayOption>::const_iterator itArrayOpt = arrayOpt.begin();
            map<string, Plot::SaveArrayOption>::const_iterator itArrayOptEnd = arrayOpt.end();
            string genErrMsg("Required tclout plot syntax:\n");
            genErrMsg += "    tclout plot <option> <array> [<plot group no.>]\n\n";
            genErrMsg += "where <option> is a valid plot option and <array> is one of:\n    ";
            while (itArrayOpt != itArrayOptEnd)
            {
               genErrMsg += itArrayOpt->first + "  ";
               ++itArrayOpt;
            }
            string plotCom = savPlotCom;
            Plot::SaveArrayOption arrayEnum = savArrayEnum;
            int plotGroup = savPlotGroup;
            bool rerouteFlag = false;
            for (size_t i=0; i<iParams.size(); ++i)
            {
               switch (iParams[i])
               {
                  case 0: // plot command
                  {
                     // Need to check for "model[:<name>]"
                     plotCom = XSparse::returnDelimitedArgument(params[i], string(":"));
                     if (plotCom.length())
                        plotCom = XSutility::lowerCase(plotCom);
                     if (params[i].length())
                        plotCom += " " + params[i];                        
                  }
                     break;
                  case 1: // array option
                  {
                     const string arrayChoice = XSutility::lowerCase(params[i]);
                     map<string, Plot::SaveArrayOption>::const_iterator itChoice = arrayOpt.lower_bound(arrayChoice);
                     if (itChoice == arrayOpt.end() || itChoice->first.find(arrayChoice) != 0)
                     {
                        tcerr << genErrMsg <<endl;
                        return TCL_ERROR;
                     }
                     arrayEnum = itChoice->second;
                     // if user enters "plot model y" they really want what we store
                     // in PlotGroup.model.  Therefore re-route it...
                     const string testForMod("model");
                     // plotCom may have more than 1 word due to named model.
                     string tmpPlotCom(plotCom);
                     const string upToSpace = XSparse::returnDelimitedArgument(tmpPlotCom, string(" "));
                     if (testForMod.find(upToSpace) == 0 && arrayEnum == Plot::Y_AR)
                     {
                        arrayEnum = Plot::MODEL_AR;
                        rerouteFlag = true;
                     }
                  }
                     break;
                  case 2: // plot group number
                  {
                     istringstream iss(params[i]);
                     if (!(iss >> plotGroup) || !iss.eof() || plotGroup < 1)
                     {
                        tcerr << "Improper plot group number specifier\n"
                           << genErrMsg << endl;
                        return TCL_ERROR;
                     }
                  }
                     break;
                  default:
                     break;
               } // end switch/case
            } // end Params loop
            pair<Plot::SaveArrayOption,int> arrayInfo(arrayEnum, plotGroup);
            XSContainer::plot->saveArrayInfo(arrayInfo);
            string tclPlotCom("plot ");
            tclPlotCom += plotCom;
            // Exceptions thrown during plot should be caught and handled
            // in xsPlot. Only an error code should exist by this point. 
            int plotStatus = Tcl_Eval(interp, tclPlotCom.c_str());
            if (plotStatus == TCL_OK)
            {
               const vector<Real>& outputArray = XSContainer::plot->savedPlotArray();
               // outputArray size will be 0 if plot was unable to find
               // what user requested.
               size_t sz = outputArray.size();
               if (sz)
               {
                  for (size_t i=0; i<sz; ++i)
                  {
                     ss << outputArray[i] << "   ";
                  }
               }
               else ss << .0;
               savPlotCom = plotCom;
               savArrayEnum = arrayEnum;
               if (rerouteFlag)  savArrayEnum = Plot::Y_AR;
               savPlotGroup = plotGroup;
            }
            else
            {
               ss << .0;
            }
         }
            break;
         case PLOTGRP: // plotgrp
            ss << datasets->numberOfPlotGroups();
            break;
         case QUERY: // query
            if (fit->queryMode() == Fit::ON)
            {
               outToTcl = "on";
            }
            else if (fit->queryMode() == Fit::YES)
            {
               outToTcl = "yes";
            }
            else
            {
               outToTcl = "no";
            }
            ss << outToTcl;
            break;
         case RATE: //rate
            // tclout rate has two options: rate or rate n, according
            // to whether the user wants the total count rate or the
            // count rate for a specific data set and model.
            if (objc < 3)
            {
                    tcerr << " syntax: tclout rate <n | all >\n";
                    return TCL_ERROR;
            }
            else
            {
               string arg = XSutility::lowerCase(string(Tcl_GetString(objv[2])));
               if (arg == "all")
               {
                  specNum = 0;
               }
               else
               {
                  if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
                  {
                     tcerr << " syntax: tclout rate <n | all >\n";
                     return TCL_ERROR;
                  }
               }
               
               // If specNum = 0, countRate function sums over all spectra
               Real modelRate = models->countRate(specNum);
               Real dataRate(.0);
               Real variance(.0);
               Real totFlux(.0);
               Real percent(.0);
               if (specNum)
               {
                  const SpectralData* sd = datasets->lookup(specNum);
                  dataRate = sd->netFlux();
                  variance = sd->netVariance();
                  totFlux = sd->totalFlux();
               }
               else  // sum over all spectra
               {
                  const size_t N = datasets->numberOfSpectra();
                  for ( size_t j = 1; j <= N; ++j)   
                  {
                          SpectralData* spectrum = datasets->lookup(j); 
                          dataRate += spectrum->netFlux();
                          totFlux += spectrum->totalFlux();
                          Real curVar = spectrum->netVariance();
                          variance +=  curVar*curVar;
                  }
                  variance = sqrt (variance);                         
               }
               if (totFlux != .0)
                  percent = 100.*dataRate/totFlux;
               ss << dataRate <<" "<< variance <<" "<< modelRate << " " << percent;;
            }
            break;
         case RESPONSE: // response n
            if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
            {
               return TCL_ERROR;               
            }
            else
            {
               SpectralData* sd = datasets->lookup(specNum);
               string respName("");
               const std::vector<Response*>& resps = sd->detector();
               for (size_t i=0; i<resps.size(); ++i)
               {
                  const RealResponse* rrsp = dynamic_cast<RealResponse*>(resps[i]);
                  if (rrsp)
                  {
                     respName += rrsp->rmfName();
                     respName += " ";
                  }
                  else
                  {
                     const MultiResponse* mrsp = dynamic_cast<MultiResponse*>(resps[i]);
                     if (mrsp)
                     {
                        const StringArray& rmfNames = mrsp->rmfNames();
                        for (size_t j=0; j<rmfNames.size(); ++j)
                        {
                           respName += rmfNames[j];
                           respName += " ";
                        }
                     }
                  }
               }
               ss << respName;
            }
            break;            
         case SIMPARS:  // simpars
            if (!fit->isStillValid())
            {
               tcerr << "An up-to-date fit is first required to randomize model parameters."<<endl;
               return TCL_ERROR;
            }
            else
            {
               const std::map<int,ModParam*>& varPars = fit->variableParameters();
               std::map<int,ModParam*>::const_iterator it = varPars.begin();
               std::map<int,ModParam*>::const_iterator itEnd = varPars.end();
               std::vector<Real> savePars;
               while (it != itEnd)
               {
                  savePars.push_back(it->second->value());
                  ++it;
               }
               bool isErr = false;
               try
               {
                  fit->reportErrorSimulationMethod();
                  // this may throw
                  fit->randomizeModelParameters(true);
                  // This should retrieve all active/on models in 
                  // alphabetical order, and internally by data 
                  // group number.
                  const std::map<string,bool>& modNames = models->activeModelNames();
                  std::map<string,bool>::const_iterator itNames = modNames.begin();
                  while (itNames != modNames.end())
                  {
                     if (itNames->second)
                     {
                        const string& modName = itNames->first;
                        const Model* firstMod = models->lookup(modName);
                        if (!firstMod->isActive()) 
                        {
                           ++itNames;
                           continue;
                        }
                        // At this point, modName must refer to 
                        // active/on models.
                        size_t nParsInMod = firstMod->numberOfParameters();
                        // This is O(N^2) lookup to get the mods in order
                        // of data group number...BAD!!!  But, as N is
                        // never likely to be > 10^3, this should be
                        // OK....for now.
                        size_t N = datasets->getNumberOfGroupsForSource(firstMod->sourceNumber());
                        for (size_t iGroup=1; iGroup<=N; ++iGroup)
                        {
                           for (size_t iPar=1; iPar<=nParsInMod; ++iPar)
                           {
                              Real parVal = models->lookupParameter(iPar+
                                           (iGroup-1)*nParsInMod,modName)->value();
                              ss << parVal << " ";
                           }
                        }
                     }
                     ++itNames;
                  }
               }
               catch (YellowAlert&)
               {
                  isErr = true;
               }
               // Restore original parameters and fit state.
               it = varPars.begin();
               size_t nVarPars = varPars.size();
               for (size_t i=0; i<nVarPars; ++i, ++it)
               {
                  it->second->setValue(savePars[i],'v');
               }
               fit->statMethod()->perform(fit);
               fit->isStillValid(true);
               if (isErr)
                  return TCL_ERROR;
            }
            break;
         case SOLAB: // solab
            {
               const string& abund = FunctionUtility::ABUND();
               const vector<float>& abundances = 
                        FunctionUtility::abundanceVectors(abund);
               const size_t& N = FunctionUtility::NELEMS();
               
               std::ostringstream oss;
               oss.precision(3);
               oss.setf(ios_base::scientific | ios_base::showpoint);
               for ( size_t j = 0; j < N; ++j)
               {
                  oss << abundances[j] << " ";
               }
               ss << oss.str();
            }
            break;
         case STAT: // stat
               ss << fit->statistic();
            break;
         case STATMETHOD: 
               ss << fit->statMethod()->fullName();
            break;
         case MARGIN:  // margin probability | [<modName>:]<parameterNum>
         case STEPPAR:  // steppar statistic | [<modName>:]<parameterNum>
         {
            const string cmdLabel = (selected == MARGIN) ? string("margin") :
                                string("steppar");
            const string optLabel = (selected == MARGIN) ? string("probability") :
                                string("statistic");
            if (objc < 3)
            {
               tcerr << " Proper syntax is: tclout " << cmdLabel <<" "
                        << optLabel << " | [<modName>:]<parameterNum>" << std::endl;
               return TCL_ERROR;
            }
            if (selected == STEPPAR && !fit->isStillValid())
            {
               tcerr << " ***Error: A valid fit does not exist.  Fit and steppar must be (re)run\n"
                    << "   before values can be sent to tclout steppar." << std::endl;
               // To prevent confusion, clear out old tclout.
               ss << " ";
               Tcl_SetVar(interp, XSPEC_TCLOUT, ss.str().c_str(), 0); 
               return TCL_ERROR;                 
            }
            
            string testOpt = XSutility::lowerCase(string(Tcl_GetString(objv[2])));
            // For either command, we only need what's in the
            // Grid base class.
            const Grid* grid = (selected == MARGIN) ? 
                        static_cast<Grid*>(fit->chainManager()->marginGrid()) :
                        static_cast<Grid*>(fit->stepGrid());
            if (!grid)
            {
               tcerr << " No " << cmdLabel << " results currently exist." << std::endl;
               return TCL_ERROR;
            }
            if (optLabel.find(testOpt) != string::npos)
            {
               // The easier case: simply return the fit statistic values.
               const RealArray& gridVals = grid->getGrid();
               size_t nVals = gridVals.size();
               if (!nVals)
               {
                  // Not sure if this case is even possible, but this can't hurt.
                  ss << "  ";
               }
               for (size_t i=0; i<nVals; ++i)
               {
                  ss << gridVals[i] << "  ";
               }
            }
            else
            {
               // iPar is not the par index number (either user entered or
               // internal).  It merely refers to its offset in the grid
               // parameter array (ie. 2nd steppar param has iPar = 1).
               size_t iPar = 0;
               size_t nIntervalOffset = 0;
               const Grid::ParameterSpec* foundPar = 0;
               Grid::SpecContainer::const_iterator itPars = grid->getParameter().begin();
               Grid::SpecContainer::const_iterator itParsEnd = grid->getParameter().end();
               if (selected == STEPPAR)
               {
                  // For steppar, parameter arrays have 1 more value than
                  // interval settings, so...
                  nIntervalOffset = 1;
                  Parameter* paramObj = HandlerUtils::getFromStringInt<Parameter>(objc, objv);
                  if (!paramObj)
                  {
                     tcerr << "***Error: Unable to retrieve parameter." << std::endl;
                     return TCL_ERROR;
                  }
                  ModParam* modPar = dynamic_cast<ModParam*>(paramObj);
                  if (!modPar)
                  {
                     tcerr << " ***Error: Specified parameter is not a model parameter." 
                          << std::endl;
                     return TCL_ERROR;
                  }
                  while (!foundPar && itPars != itParsEnd)
                  {
                    if ((*itPars)->address == modPar)
                    {
                       foundPar = *itPars;
                    }
                     ++iPar;
                     ++itPars;
                  }
               }
               else
               {
                  string modName;
                  size_t parNum = string::npos;
                  if (HandlerUtils::verifyStringInt<Parameter>(objc, objv,
                                3, modName, parNum))
                  {
                     while (!foundPar && itPars != itParsEnd)
                     {
                        if ((*itPars)->name == modName &&
                              static_cast<size_t>((*itPars)->index) == parNum)
                        {
                           foundPar = *itPars;
                        }
                        ++iPar;
                        ++itPars;
                     }
                  }
               }
               if (foundPar)
               {
                  // If we're here iPar must be >= 1, so decrement on a 
                  // size_t is OK.
                  --iPar;
                  const std::vector<Real>& parVals = foundPar->parameterValues;
                  const size_t nParVals = parVals.size();
                  const size_t nGrid = grid->getGrid().size();
                  std::vector<Real> repetitiveVals(nGrid, .0);
                  size_t innerRepCount = 1;
                  for (size_t i=0; i<iPar; ++i)
                  {
                     innerRepCount *= (grid->getParameter()[i])->intervals 
                                        + nIntervalOffset;
                  }
                  size_t outerRepCount = nGrid/(innerRepCount*nParVals);
                  for (size_t i=0; i<outerRepCount; ++i)
                  {
                     size_t iOffset = i*nParVals*innerRepCount;
                     for (size_t j=0; j<nParVals; ++j)
                     {
                        Real parVal = parVals[j];
                        size_t jOffset = j*innerRepCount;
                        for (size_t k=0; k<innerRepCount; ++k)
                        {
                           repetitiveVals[k + iOffset + jOffset] = parVal;
                        }
                     }
                  }
                  // Success
                  for (size_t i=0; i<nGrid; ++i)
                  {
                     ss << repetitiveVals[i] << "  ";
                  }                                         
               }
               else
               {
                  tcerr << " ***Error: Parameter " << string(Tcl_GetString(objv[2]))
                      << " not found in results of last " << cmdLabel << "." 
                      << std::endl;
                  return TCL_ERROR;
               }
            } // end handling [<modName>:]<parNum>
         }          
            break;        
         case VARPAR: // varpar
            if (fit->isStillValid())
            {
               size_t nVar = fit->variableParameters().size();
               ss << nVar;
            }
            else
            {
               tcerr << " Valid fit does not exist." << std::endl;
               return TCL_ERROR;
            }
            break;
         case WEIGHT:
            ss << DataUtility::statWeight().name();
            break;
         case XFLT: // xflt n
            if (!HandlerUtils::getSpecNumParameter(objc, objv, opt, 3, specNum))
            {
               return TCL_ERROR;               
            }
            else
            {
               const SpectralData* sd = datasets->lookup(specNum);
               const std::map<size_t, Real>& vals = sd->xflt();
               size_t sz = vals.size();
               ss << sz;
	       
	       std::map<size_t, Real>::const_iterator
		   i_xfltBeg = vals.begin(),
		   i_xfltEnd = vals.end();

	       while(i_xfltBeg != i_xfltEnd)
	       {
		   ss << " " << i_xfltBeg->second;
		   ++i_xfltBeg;
	       }

            }
            break;
         default:
            {
               string msg("TCLOUT option ");
               msg += opt;
               msg += " not yet implemented";
               tcout << msg << std::endl; 
               XSutility::printValidOptions(tcout, cmd, option);
            }
            break;
      } // end switch
      Tcl_SetVar(interp, XSPEC_TCLOUT, ss.str().c_str(), 0);                  
   }
   else
   {
      string errMsg(inputOpt);
      errMsg += " is not currently a valid TCLOUT option.\n";
      tcerr << errMsg << std::endl;
      XSutility::printValidOptions(tcout, cmd,option);
      return TCL_ERROR;
   }
   
   return TCL_OK;
}