//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: nupower.cc,v 1.9 2004-09-10 09:54:40 cmv Exp $
//
// Description:
//      Aim of the class: To give the energy density for a nupower
//                        The unity used here is W/m^2/Hz/sr
//
// History (add to end):
//      Sophie   Oct, 1999  - creation
//
//------------------------------------------------------------------------

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

#include "nupower.h"


/*! 
 * \class SOPHYA::PowerLawSpectra 
 \ingroup SkyT
 * This class corresponds to a power law radiation spectrum.
*/
PowerLawSpectra::PowerLawSpectra()
{
}
/*! Constructor: The arguments corresponds to the Power Law Spectrum
  equation:
<a name="psldef"> </a>
\f[
\hbox{flux}(\nu) = A ({\nu-\nu_0\over \delta\nu})^B
\f]
*/
PowerLawSpectra::PowerLawSpectra(double A, double B, double nu0, double dnu, double numin, double numax)
  : RadSpectra(numin, numax)
{
  _a = A;
  _b = B;
  _nu0 = nu0;
  _dnu = (dnu > 1.e-19) ? dnu : 1.;
}


PowerLawSpectra::~PowerLawSpectra()
{
}
/*! Flux function according to the <a href="#psldef"> above </a>
 equation */
double
PowerLawSpectra::flux(double nu) const
{
  if ((nu< _numin) || (nu > _numax)) return(0.);
  double tmp = (nu-_nu0)/_dnu;
  if (tmp <= 0.)  return(0.); 
  else return(  _a * pow(tmp, _b) );
}



void
PowerLawSpectra::Print(ostream& os) const
{
  os << "PowerLawSpectra::Print f = a ((nu-nu0)/dnu)^b " << _a << "((nu-" << _nu0 
     << ") / " << _dnu << ") ^ " << _b << endl;
  os << " - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
  os << "MeanFreq= " << meanFreq() << "  Emission= " << flux(meanFreq()) << endl;
  os << "PeakFreq= " << peakFreq() << "  Emission= " << flux(peakFreq()) << endl;

}


/*
void 
PowerLawSpectra::WriteSelf(POutPersist& s)  
{
  s.PutR8(this->giveNorm());
  s.PutR8(this->giveNu0());
  s.PutR8(this->giveDNu());
  s.PutR8(this->giveExp());
  s.PutR8(this->minFreq());
  s.PutR8(this->maxFreq());
}

void
PowerLawSpectra::ReadSelf(PInPersist& s)  
{
  s.GetR8(_a);
  s.GetR8(_nu0);
  s.GetR8(_dnu);
  s.GetR8(_b);
  s.GetR8(_numin);
  s.GetR8(_numax);
  cout << " Norm - Nu0 - DNu - Exp - minFreq - maxFreq " << endl; 
  cout << _a << "-" << _nu0  << "-" << _dnu << "-" << _b << "-" << _numin << "-" << _numax << endl;
}

*/