//   Read the documentation to learn more about C++ code generator
//   versioning.
//	  %X% %Q% %Z% %W%

// Plot
#include <XSPlot/Plot/Plot.h>
// DataContainer
#include <XSModel/GlobalContainer/DataContainer.h>
// ModelContainer
#include <XSModel/GlobalContainer/ModelContainer.h>
// PlotUnfolded
#include <XSPlot/Commands/PlotUnfolded.h>
#include <XSPlot/Plot/PlotGroup.h>
#include <XSUtil/Numerics/Numerics.h>
#include <XSContainer.h>
#include <XSsymbol.h>


// Class PlotUnfolded 

PlotUnfolded::PlotUnfolded (Plot* plot, const std::string& name)
  : PlotCommand(plot,name)
{
  modelRequired(Plot::INACTIVE);
}


void PlotUnfolded::primaryGraph (const string& value)
{
  Plot* pPlot = plot();

  bool constructAddComponents(true);

  Plot::XaxisMode saveXoption (pPlot->xOption());
  if ( saveXoption == Plot::CHANNELS )
  {
     pPlot->setXOption("energy");
  }
  // Always plot log scale for x and y axis for model related plots.
  pPlot->plotXLog(true);
  pPlot->plotYLog(true);

  try
  {
     pPlot->setSelectedLabel(2,"tu");
     if ( name() == "ufspec")
     {
           pPlot->autoSelectedLabels("u*");
     }
     else if ( name() == "eufspec")
     {
           pPlot->autoSelectedLabels("v*");       
     }
     else if ( name() == "eeufspec")
     {
           pPlot->autoSelectedLabels("w*");       
     }
     if ( !pPlot->showAddComponents())
     {
           // if this hasn't already been done, do it here.       
  	   XSContainer::models->clearSources();
  	   XSContainer::models->makeSourceComponents();
     }      

     XSContainer::models->foldSources();

     pPlot->makeBinnedPlotArrays(true);

     fillDataAttributes(pPlot->dataStyle(),0,2);
     fillModelAttributes();
     fillSourceAttributes();
     manipulate();

     pPlot->fixYMax(pPlot->ranges("y").second*2.);
     pPlot->fixYMin(std::max(pPlot->ranges("y").first,pPlot->ranges("y").second*1.E-5));
     pPlot->rangeDefaultFlags(15);


   // graph it.

     pPlot->bundlePlotVectors(constructAddComponents);
     pPlot->prepareTheGraph();
     pPlot->makeTheGraph(); 
  }
  catch (YellowAlert &)
  {
     if ( saveXoption == Plot::CHANNELS ) pPlot->setXOption("channel");
     throw;
  }
  if ( saveXoption == Plot::CHANNELS ) pPlot->setXOption("channel");
}

void PlotUnfolded::manipulate ()
{
  PlotGroupContainer::iterator  g = plot()->groups().begin();
  PlotGroupContainer::iterator  gEnd = plot()->groups().end();

  typedef std::vector<Real> RealVector;

  const string Y("y");

  plot()->ranges(Y).first = LARGE;
  plot()->ranges(Y).second = -LARGE;

  bool modelStep(true);
  bool sourceStep(true);
  setModelLineStep(modelStep,sourceStep);

  while (g != gEnd)
  {
     PlotGroup* gr = *g;
     int arraySize = gr->n;
     const size_t nMods = gr->model.size();
     RealVector& X      = gr->xAxis.data;
     RealVector& dX      = gr->xAxis.errors[0];
     RealVector& yData  = gr->yData.data;
     RealVector& yError = gr->yData.errors[0];
     RealVector& foldedModel  = gr->model[0].data;
     RealVector& unfoldedModel = gr->auxData[0].data;

     setDataColor(gr,PlotStyle::Colour(gr->plotGroup)); 
     setModelColor(gr,PlotStyle::Colour(gr->plotGroup));  
     setSourceColor(gr,PlotStyle::Colour(gr->plotGroup)); 
     if (nMods > 1)
     {
        // This will contain sum of models, unfolded.
        unfoldedModel.resize(arraySize,.0);
        for (int l = 1; l < nMods; ++l)
        {
           RealVector& currentUnfolded = gr->auxData[l].data;
           for ( int k = 0; k < arraySize; ++k)
           {
                   unfoldedModel[k] += currentUnfolded[k];       
           }
          // setModelColor(gr,PlotStyle::Colour(gr->plotGroup+2 + l),l);  
         //  setSourceColor(gr,PlotStyle::Colour(gr->plotGroup+2 + l),l); 
        }

     }
     for ( int k = 0; k < arraySize; ++k)
     { 
        // Apply uf/f ratio to data.
        if ( foldedModel[k] != 0.0)
        {                        
              Real multiplier (unfoldedModel[k]/foldedModel[k]);
              yData[k]  *= multiplier;
              yError[k] *= multiplier;
        }
        else
        {
              yData[k]  = Plot::NODATA();
              yError[k] = Plot::NODATA();                    
        }
     }

     // Move all unfolded arrays from their workspaces into what up to now
     // was holding folded.
     for (int i=0; i<gr->model.size(); ++i)
     {
        gr->model[i].data = gr->auxData[i].data;
     }
     for (int i=0; i<gr->sources.size(); ++i)
     {
        PlotVectorList::iterator ucs = gr->sourceFluxes[i].begin();
        PlotVectorList::iterator ucsEnd = gr->sourceFluxes[i].end();
        PlotVectorList::iterator cs = gr->sources[i].begin();
        while (ucs != ucsEnd)
        {
           cs->data = ucs->data;
           ++ucs, ++cs;
        }
     }

     RealVector multiplier(arraySize,1.);
     if (name() == "eufspec")
     {
        if (plot()->xOption() == Plot::WAVELENGTH)
        {
           // For wavelength case, still want to multiply by E, though
           // X and dX are lengths (Angstroms).
           for (int k=0; k<arraySize; ++k)
           {
              Real commonFact = Numerics::KEVTOA/(X[k]*X[k]+dX[k]*dX[k]);
              Real E = X[k]*commonFact;
              Real dE = dX[k]*commonFact;
              if (plot()->waveUnitsHz())
                 multiplier[k] = Numerics::KEVTOJY*sqrt((E - dE)*(E + dE));
              else
                  multiplier[k] = sqrt((E - dE)*(E + dE));  
           }
        }
        else
        {
           for (int k=0; k<arraySize; ++k)
           {
              multiplier[k] = sqrt((X[k] - dX[k])*(X[k] + dX[k]));
           }
        }
     }
     else if (name() == "eeufspec")
     {
        if (plot()->xOption() == Plot::WAVELENGTH)
        {
           for (int k=0; k<arraySize; ++k)
           {
              Real commonFact = Numerics::KEVTOA/(X[k]*X[k]+dX[k]*dX[k]);
              Real E = X[k]*commonFact;
              Real dE = dX[k]*commonFact;
              if (plot()->waveUnitsHz())
                 multiplier[k] = Numerics::KEVTOERG*Numerics::KEVTOHZ*(E - dE)*(E + dE); 
              else
                 multiplier[k] = (E - dE)*(E + dE);  
           }
        }
        else
        {
           for (int k=0; k<arraySize; ++k)
           {
              multiplier[k] = (X[k] - dX[k])*(X[k] + dX[k]);
           }
        }
     }

     for ( int k = 0; k < arraySize; ++k)
     { 
        if ( multiplier[k] != 1.)
        {
           Real mult = multiplier[k];       
           yData[k] *= mult;
           yError[k] *= mult;

           for (int i=0; i<gr->model.size(); ++i)
           {
              gr->model[i].data[k] *= mult;
           }
           for (int i=0; i<gr->sources.size(); ++i)
           {
              PlotVectorList::iterator cs = gr->sources[i].begin();
              PlotVectorList::iterator csEnd = gr->sources[i].end();
              while (cs != csEnd)
              {
                 cs->data[k] *= mult;
                 ++cs;
              }
           }
        } // end if multiplier[k] != 1.
     }  // end k arraySize loop

     // set range values

     for ( int k = 0; k < arraySize; ++k)
     { 
       if ( yData[k] != Plot::NODATA() ) {
           plot()->ranges(Y).first  
                     = std::min(yData[k],plot()->ranges(Y).first);
           plot()->ranges(Y).second  
                     = std::max(yData[k],plot()->ranges(Y).second);
       }
     }


     ++g;           
  } // end plot group loop
}

// Additional Declarations