source: Sophya/trunk/SophyaLib/SkyT/nupower.cc@ 668

//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: nupower.cc,v 1.4 1999-11-29 14:16:05 ansari 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 "machdefs.h"
#include <iostream.h>
#include <math.h>

#include "nupower.h"

//----------------
// Constructor --
//----------------
PowerLawSpectra::PowerLawSpectra()
{
}

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()
{
}

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::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;
}


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;

}