source: Sophya/trunk/SophyaLib/SkyT/blackbody.cc@ 3865

//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: blackbody.cc,v 1.9 2004-09-10 09:54:39 cmv Exp $
//
// Description:
//      Aim of the class: To give the energy density for a blackbody
//                        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 "blackbody.h"

/*! 
 * \class SOPHYA::BlackBody 
 \ingroup SkyT
 This class corresponds to the emission spectrum of a
 blackbody radiation.
 */
//----------------
// Constructor --
//----------------
/*! Constructor: needs a temperature. Otherwise set to ConvTools::tcmb */
BlackBody::BlackBody(double temperature)
  : RadSpectra(10., 10000.)
{
  _temperature = temperature;
}


BlackBody::~BlackBody()
{
}
/*! Black Body Flux Function:
  \f[
  I_\nu = {2 h_{pl} (1.10^9*\nu)^3 \over c^2 (e^{{h_{pl}(1.10^9*\nu) \over kT}} -1)}
  \f]
*/
double
BlackBody::flux(double nu) const
{
  if(nu < -1.e99) nu = -1.e99;
  if(nu > 1.e99) nu = 1.e99;
  double temperature = getTemperature();
  if(nu==0.) return 0.;
  double hpl = ConvTools::hpl;
  double cel = ConvTools::cel;
  double kb  = ConvTools::kb;
  double puiss1 = nu*pow(10.,9);
  if(puiss1 >  1.e99)   puiss1=1.e99;
  if(puiss1 < -1.e99)   puiss1=-1.e99;
  double puiss2 = hpl*nu*pow(10.,9)/(kb*temperature);
  if(puiss2 >  1.e99)   puiss2=1.e99;
  if(puiss2 < -1.e99)   puiss2=-1.e99;
  
  double result=
   (2*hpl* pow( puiss1 ,3))/( pow(cel,2)*(exp(puiss2)-1));
  return result;
}


void
BlackBody::Print(ostream& os) const
{
  os << "BlackBody::Print Temp= " << getTemperature() 
     << " - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
  os << "MeanFreq= " << meanFreq() << "  Emission= " << flux(meanFreq()) << endl;
  os << "PeakFreq= " << peakFreq() << "  Emission= " << flux(peakFreq()) << endl;

}

/*
void 
BlackBody::WriteSelf(POutPersist& s)  
{
  s.PutR8(this->getTemperature());
  s.PutR8(this->minFreq());
  s.PutR8(this->maxFreq());
}

void
BlackBody::ReadSelf(PInPersist& s)  
{
  s.GetR8(_temperature);
  s.GetR8(_numin);
  s.GetR8(_numax);
  cout << " Temperature - minFreq - maxFreq " << endl; 
  cout << _temperature << "-" << _numin << "-" << _numax << endl;
}

*/