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

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


// Class PlotModel 

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


void PlotModel::manipulate ()
{
  const string Y("y");

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

  PlotGroupContainer::iterator  g = plot()->groups().begin();
  PlotGroupContainer::iterator  gEnd = plot()->groups().end();
   while (g != gEnd)
   {
      PlotGroup* gr = *g;
      int arraySize = gr->n;
      std::vector<Real>& X      = gr->xAxis.data;
      std::vector<Real>& dX      = gr->xAxis.errors[0];
      std::vector<Real> multiplier(arraySize);
      if ( name() == "model" )
      {
	if (plot()->xOption() == Plot::WAVELENGTH && plot()->waveUnitsHz() == 1)
	{
	    for ( int k = 0; k < arraySize; ++k)
	    {
	        Real dHz = Numerics::KEVTOHZ*Numerics::KEVTOA*dX[k]/(X[k]*X[k]-dX[k]*dX[k]);
	      	multiplier[k] = 1.0/(2*dHz);
	    }         
	} else {
	    for ( int k = 0; k < arraySize; ++k)
	    {
		multiplier[k] = 1./(2*dX[k]);       
	    }         

	}                   
      }
      else
      {         
         if ( name() == "emodel" )
         {
            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() == 1)
                  {
		     Real dHz = Numerics::KEVTOHZ*Numerics::KEVTOA*dX[k]/(X[k]*X[k]-dX[k]*dX[k]);
                     multiplier[k] = Numerics::KEVTOJY*sqrt((E - dE)*(E + dE))/(2*dHz); 
                  }
                  else
                     multiplier[k] = sqrt((E - dE)*(E + dE))/(2*dX[k]);  
               }
            }
            else
            {
               for ( int k = 0; k < arraySize; ++k)
               {
                  multiplier[k] = sqrt((X[k] - dX[k])*(X[k] + dX[k]))/(2*dX[k]);   
               }
            }                                                    
         }
         else if ( name() == "eemodel" )
         {
            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() == 1)
                  {
		     Real dHz = Numerics::KEVTOHZ*Numerics::KEVTOA*dX[k]/(X[k]*X[k]-dX[k]*dX[k]);
                     multiplier[k] = Numerics::KEVTOERG*Numerics::KEVTOHZ*(E - dE)*(E + dE)/(2*dHz); 
                  }
                  else
                     multiplier[k] = (E - dE)*(E + dE)/(2*dX[k]);
               }
            }
            else
            {
               for ( int k = 0; k < arraySize; ++k)
               {
                  multiplier[k] = (X[k] - dX[k])*(X[k] + dX[k])/(2*dX[k]);   
               }
            }                                                    
         }
      }
      for (size_t j = 0; j < gr->model.size(); ++j)
      {
         std::vector<Real>& model  = (*g)->model[j].data;
         std::vector<Real>* modelErr = 0;
         if ((*g)->model[j].errors.size())
         {
            modelErr = &(*g)->model[j].errors[0];
         }
         bool error = false;
         if (modelErr && modelErr->size() == model.size())
         {
            error = true;
         }
         for ( int k = 0; k < arraySize; ++k)
         {
            model[k] *= multiplier[k];  
            plot()->ranges(Y).first 
                    = std::min(model[k],plot()->ranges(Y).first);
            plot()->ranges(Y).second 
                    = std::max(model[k],plot()->ranges(Y).second);
            if (error) (*modelErr)[k] *= multiplier[k];     
         }
      }
      for (size_t j = 0; j < gr->sources.size(); ++j)
      {
         PlotVectorList::iterator sl = gr->sources[j].begin();
         PlotVectorList::iterator slEnd = gr->sources[j].end();
         while ( sl != slEnd )
         {
            PlotVector& s = *sl;
            for ( int k = 0; k < arraySize; ++k) s.data[k] *= multiplier[k];
            ++sl;       
         }
      }
      ++g;       
   }          
}

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

  std::vector<Model*> mods = XSContainer::models->lookupModelGroup(value);
  if (!mods.size())
  {
     string msg;
     if (value == Model::DEFAULT())
     {
        msg = "No model has been entered with the default name";
     }
     else
     {
        msg = "No model has been entered with name = " + value;
     }
     throw Plot::ModelNotFound(msg);
  }

// manufacture necessary data arrays into pPlot's PlotGroup container.
// For this particular graph need X axis set by setXaxis using the 
// Response's energy array (Ehi, Elo, not EboundsMax, EboundsMin)
// and its efficiency array.

  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
  {  
     // set titles and labels.

     pPlot->setSelectedLabel(2,"tm");
     if ( name() == "model")
     {
           pPlot->autoSelectedLabels("m*");
     }
     else if ( name() == "emodel")
     {
           pPlot->autoSelectedLabels("n*");       
     }
     else if ( name() == "eemodel")
     {
           pPlot->autoSelectedLabels("o*");       
     }

     // if this hasn't already been done, do it here.       
     XSContainer::models->clearSources();
     //XSContainer::models->makeSourceComponents();

     // Note: ModelContainer's makeSourceComponents, which
     // is already called in commonInit, does not
     // currently perform on inactive models.  Therefore, it
     // is called below singly on each model.
     for(int i = 0; i < mods.size(); ++i)
	 mods[i]->makeSourceComponents();

     pPlot->makePlotArraysFromModels(value, true);


     manipulate();
     fillModelAttributes(1,2);
     fillSourceAttributes();
     // Default range will only show at most 3 orders of magnitude.  These
     // adjustments assume ymin,ymax are positive, but then we're already 
     // making that assumption by forcing ylog on.
     Real ymax = pPlot->ranges("y").second;
     Real ymin = std::max(pPlot->ranges("y").first,ymax*1.E-3);
     // Add an epsilon to range to deal with case where ymax = ymin.
     const Real epsRange = .01; 
     pPlot->fixYMin(ymin*(1.0 - epsRange));
     pPlot->fixYMax(ymax*(1.0 + epsRange));
     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");
}

// Additional Declarations