source: Sophya/trunk/SophyaLib/SkyT/radspec.cc@ 3613

//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: radspec.cc,v 1.8 2004-09-10 09:54:40 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 "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <typeinfo>
#include <math.h>

#include "radspec.h"
#include "integ.h"
/*!
   \defgroup SkyT SkyT module
   This module contains classes and functions which define 
   several radiation spectra and filter responses
*/

/*! 
 * \class SOPHYA::RadSpectra
 * \ingroup SkyT
 * This class is an abstract base class for radiation emission spectra. The flux() function returns the value of the flux (the spectral      <BR>
 * energy distribution) as a function of the frequency.  As in the SpectralResponse class, the () operator has been redefined <BR>
 * at this level, so that the user can access the flux value, either by  calling the function or directly by using this operator. <BR>
 * For all the sub-classes, \nu is given in units of Hz and
 * the flux is returned in units of W/m^2/sr/Hz.
*/


//----------------
// Constructor --
//----------------
/*! Default constructor */
/*!
  The constructor takes as an argument the minimum
  and the maximum frequency of the spectrum, if any. <BR>
  In the case the user does not want to specify these
  values, there are set respectively to 0. and 9.E49
  by default.
*/
RadSpectra::RadSpectra(double numin, double numax)
{
  _numin = numin;
  _numax = numax;
}


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

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

double 
RadSpectra::minFreq()  const
{
  return _numin;
}

double 
RadSpectra::maxFreq()  const
{
  return _numax;
}

double 
RadSpectra::meanFreq()  const
{
  double result = (_numax+_numin)/2.;
  return result;
}


/* 
   The peakFreq() function returns the value of the 
   frequency for the maximum value of the flux
*/
double 
RadSpectra::peakFreq()  const
{
  double maxAnswer = -1.e99;
  double maxNu = -10;
  double nu;
  for (int i=1; i<1000;i++)
    {
      nu=(_numax-_numin)*i/1000.+_numin;
      double lookForMax = flux(nu);
      if(maxAnswer <= lookForMax) {
        maxAnswer= lookForMax;
        maxNu    = nu;
      }
    }
  return maxNu;  
}

void
RadSpectra::setMinMaxFreq(double numin, double numax)
{
  _numin = numin;
  _numax = numax;
}

// the RadSpectra_fluxFunction function is used to call TrpzInteg double(double) 
// (integration over a range of frequencies)
static RadSpectra* _raypourfinteg = NULL;
static double RadSpectra_fluxFunction(double nu)
{
   return(_raypourfinteg->flux(nu));   
}
/*! 
  The integratedFlux() function performs the integration
  of the flux function in a frequency range <BR> defined by
  f1 and f2.
*/
double 
RadSpectra::integratedFlux(double f1, double f2)  const
{
  if(f1 < this->minFreq()) f1 = this->minFreq();
  if(f2 > this->maxFreq()) f2 = this->maxFreq();
  _raypourfinteg = const_cast<RadSpectra *>(this);
  TrpzInteg I(RadSpectra_fluxFunction , f1, f2);  
  double val = (double)I;
  _raypourfinteg = NULL;  // On ne peut pas faire ca avant la destruction de I
  return(val);        
}

/*!
  Same than integratedFlux() over the frequency range 
  of definition of the flux function
 */
double 
RadSpectra::integratedFlux()  const
{
  return integratedFlux(this->minFreq(),this->maxFreq());
}

// integration using the logarithm !!
// Carefull!! Base 10....
static RadSpectra* _rayIntLog = NULL;

static double RadSpectra_logFluxFunction(double tau)
{
  double value = _rayIntLog->flux(pow(10.,tau))*pow(10.,tau);
  return(value);   
}

/*! 
  The logIntegratedFlux() function performs the integration
  of the flux function in a frequency range <BR> defined by
  f1 and f2. The integration is here performed
  on the logarithm of the flux function. 
*/
double 
RadSpectra::logIntegratedFlux(double f1, double f2)  const
{
  if(f1 < this->minFreq()) f1 = this->minFreq();
  if(f2 > this->maxFreq()) f2 = this->maxFreq();

  double f1Log = log10(f1);
  double f2Log = log10(f2);
  if(f1Log < -1.e99) f1Log = -1.e99;
  if(f2Log > 1.e99)  f2Log = 1.e99;
  _rayIntLog = const_cast<RadSpectra *>(this);
  TrpzInteg I(RadSpectra_logFluxFunction,f1Log,f2Log);
  double value = (double)I * log(10.);
  _rayIntLog = NULL;
  return(value);        
}

/*!
  same than logIntegratedFlux over the frequency range
  of definition of the flux function
 */
double 
RadSpectra::logIntegratedFlux()  const
{
  return logIntegratedFlux(_numin,_numax);
}

// the RadSpectra_filteredFlux function is used to call TrpzInteg double(double) 
// (integration over a range of frequencies with a filter)
static SpectralResponse* _filter = NULL ;
static double RadSpectra_filteredFlux(double nu) 
{
  double flux = _raypourfinteg->flux(nu);
  return(flux * _filter->transmission(nu));
}

/*! 
  The filteredIntegratedFlux() function performs the integration
  of the flux function in a frequency range <BR> defined by
  f1 and f2 convolved by a SpectralResponse filter.
*/
double 
RadSpectra::filteredIntegratedFlux(SpectralResponse const& filter, double f1, double f2) const
{
  _raypourfinteg = const_cast<RadSpectra *>(this);
  _filter = const_cast<SpectralResponse *>(&filter);
   if(f1 < this->minFreq()) f1 = this->minFreq();
   if(f2 > this->maxFreq()) f2 = this->maxFreq();

  TrpzInteg I(RadSpectra_filteredFlux,f1,f2);
  double val = (double)I;
  _raypourfinteg = NULL;
  _filter = NULL;
  return(val);        
}

/*!
  Same than filteredIntegratedFlux() over the frequency range
  defined as: <BR>
   min_freq = MAX(minfreq_flux, minfreq_filter), <BR>
   max_freq = MIN(maxfreq_flux, maxfreq_filter), <BR>
  where:
<UL>  
<LI>   minfreq_flux   is the minimum frequency of the flux definition
<LI>   maxfreq_flux   is the maximum frequency of the flux definition
<LI>   minfreq_filter is the minimum frequency of the filter definition
<LI>   maxfreq_filter is the maximum frequency of the filter definition
  </UL>
 */
double 
RadSpectra::filteredIntegratedFlux(SpectralResponse const& filter)
{
  double minOfMin = filter.minFreq();
  double maxOfMax = filter.maxFreq();
  if(minOfMin < this->minFreq()) minOfMin = this->minFreq();
  if(maxOfMax > this->maxFreq()) maxOfMax = this->maxFreq();  
  return(filteredIntegratedFlux(filter, minOfMin, maxOfMax ) );
}


// the RadSpectraVec_filteredFlux function is used to call TrpzInteg double(double) 
// (integration over a range of frequencies with a filter)
static double RadSpectra_logFilteredFlux(double tau) 
{
  double nu = pow(10.,tau);
  double flux = _raypourfinteg->flux(nu)*nu;
  double result = flux * _filter->transmission(nu);
  return(result);
}


/*! 
 *  The filteredIntegratedFlux() function performs the integration
 * of the flux function in a frequency range <BR> defined by
 * f1 and f2 convolved by a SpectralResponse filter (using the
 * logarithm of the function).
 */
double 
RadSpectra::filteredLogIntFlux(SpectralResponse const& filter, double f1, double f2) const
{
  
   _raypourfinteg = NULL;
   _filter = NULL;
   if(f1 < this->minFreq()) f1 = this->minFreq();
   if(f2 > this->maxFreq()) f2 = this->maxFreq();
   
   double f1Log = log10(f1);
   double f2Log = log10(f2);
   if(f1Log < -1.e99) f1Log = -1.e99;
   if(f2Log > 1.e99)  f2Log = 1.e99;
   _raypourfinteg = const_cast<RadSpectra *>(this);
   _filter = const_cast<SpectralResponse *>(&filter);
   TrpzInteg I(RadSpectra_logFilteredFlux,f1Log,f2Log);  
   double val = (double)I;
   _raypourfinteg = NULL;
   _filter = NULL;
   return(val* log(10.));        
}

double 
RadSpectra::filteredLogIntFlux(SpectralResponse const& filter)
{
  return(filteredLogIntFlux(filter, filter.minFreq(), filter.maxFreq() ) );
}




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

}