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

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

#include "radspec.h"
#include "integ.h"

//----------------
// Constructor --
//----------------
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;
}



// peakFreq returns the value of the frequency for the
// peak of the spectrum. 
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;  
}

// To change min-max frequency
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));   
}

double 
RadSpectra::integratedFlux(double f1, double f2)  const
{
//   cout << endl;
//   cout <<  this->minFreq() << " = " << this->maxFreq() << endl;
//   cout <<  f1 << " = " << f2 << endl;
  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);        
}
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);   
}

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

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


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

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


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;

}