//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: derivblackbody.cc,v 1.7 2004-09-10 09:54:40 cmv Exp $
//
// Description:
//      Aim of the class: To give the derivative spectrum
//                        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 "derivblackbody.h"



/*! \class SOPHYA::DerivBlackBody 
  \ingroup SkyT
 * This class corresponds to the emission spectrum of a
 * dipole (since its emission spectrum is the derivation
 * of a blackbody spectrum wrt the temperature).
 */

/*! Constructor: needs a temperature. Otherwise set to ConvTools::tcmb */
DerivBlackBody::DerivBlackBody(double temperature)
  : RadSpectra(10., 10000.)
{
  _temperature = temperature;
}


DerivBlackBody::~DerivBlackBody()
{
}

/*!
  The flux function is the derivation of the BlackBody
  flux function wrt the temperature (used e.g. for a Dipole)
  \f[
  I_\nu = {2 h_{pl} (1.10^9*\nu)^3  {h_{pl}1.10^9*\nu \over k T^2}
        {e^{{h_{pl}(1.10^9*\nu) \over kT}}\over c^2 (e^{{h_{pl}(1.10^9*\nu) \over kT}} -1)^2}}
  \f]
*/
double
DerivBlackBody::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))*(hpl*puiss1/kb)*(1/(temperature*temperature))
    *exp(puiss2)
    /(pow(cel,2)*pow(( (exp(puiss2)-1)),2));
  //  result = 1500e3*result/400e6;
  return result;
}


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

}

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

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

*/