//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: radspecvector.cc,v 1.8 2003-02-11 15:31:07 cmv Exp $
//
// Description:
//      Aim of the class: To give the energy density
//                        The unity used here is W/m^2/Hz/sr
//
// History (add to end):
//      Sophie   Oct, 1999  - creation
//
//------------------------------------------------------------------------

//---------------
// C++ Headers --
//---------------
#include "machdefs.h"
#include <iostream>
#include <math.h>
// #include <typeinfo>

#include "radspecvector.h"
#include "pexceptions.h"

//----------------
// Constructor --
//----------------

/*! 
 * \class SOPHYA::RadSpectraVec 
 \ingroup SkyT
 * One may define the radiation
 * spectrum with two vectors: <BR> one for the frequencies and the second for the
 * value of the flux function. <BR>
 * In that case, the class to use is RadSpectraVec !
 */

RadSpectraVec::RadSpectraVec()
{
}
/*! Constructor:
  \param nu is the frequency vector <BR>
  \param fdenu is the corresponding flux values (stored in a vector as well)
 */
RadSpectraVec::RadSpectraVec(Vector const & nu, Vector const & fdenu, double numin, double numax)
        : RadSpectra(numin, numax)
{
  if(nu.NElts() != fdenu.NElts())  
    throw SzMismatchError("RadSpectraVec::RadSpectraVec() - Non equal vector sizes");
  _vecOfNu = nu;
  _vecOfFDeNu = fdenu;
  _numin = nu(0);
  _numax = nu(nu.NElts()-1);
}


//--------------
// Destructor --
//--------------
RadSpectraVec::~RadSpectraVec()
{
}

//              ---------------------------
//              --  Function Definitions --
//              ---------------------------



/*! The flux function extrapolates the flux values for the 
  frequencies that are not present in the nu vector <BR>
  given at the instanciation of the class.
 */
double 
RadSpectraVec::flux(double nu) const
{
  if ( (nu < _numin) || (nu > _numax) ) return(0.);
  double value = 0.;
  int sizeVecOfNu = _vecOfNu.NElts();
  if(nu <=  _vecOfNu(0)) return _vecOfFDeNu(0);
  if(nu >=  _vecOfNu(sizeVecOfNu-1)) return _vecOfFDeNu(sizeVecOfNu-1);
  
  for (int i=1; i<sizeVecOfNu; i++)
    {
      if(nu < _vecOfNu(i))
	{
	  double up = _vecOfFDeNu(i) ;
	  double down = _vecOfFDeNu(i-1);
	  double xmin = _vecOfNu(i-1);
	  double xmax = _vecOfNu(i);
	  double a = ((up-down)/(xmax-xmin));
	  value = a*nu+(up-a*xmax);
	  return value;
	}
    }
  return _vecOfFDeNu(sizeVecOfNu-1);
}



void
RadSpectraVec::Print(ostream& os) const
{
  //  os << "RadSpectraVec::Print (" << typeid(*this).name() 
  //     << ") - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
  os << "RadSpectraVec::Print  - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
  os << "MeanFreq= " << meanFreq() << "  Emission= " << flux(meanFreq()) << endl;
  os << "PeakFreq= " << peakFreq() << "  Emission= " << flux(peakFreq()) << endl;

}

/*
void
RadSpectraVec::WriteSelf(POutPersist& s) 
{
  s.PutI4(this->NbElts());
  for (int i=0; i< this->NbElts(); i++)
    {
      s.PutR8(this->getNuVec(i));
      s.PutR8(this->getFNuVec(i));
    }
  s.PutR8(this->minFreq());
  s.PutR8(this->maxFreq());
 
}

void
RadSpectraVec::ReadSelf(PInPersist& s) 
{
  s.GetI4(_size);

  _vecOfNu.ReSize(_size);
  _vecOfFDeNu.ReSize(_size);
  for (int i=0; i< _size; i++)
    {
      s.GetR8(_vecOfNu(i));
      s.GetR8(_vecOfFDeNu(i));
    } 
  s.GetR8(_numin);
  s.GetR8(_numax);
  cout << "minFreq - maxFreq"<< endl; 
  cout <<    _numin<< "-" <<  _numax << endl;
}
*/