//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: specresp.cc,v 1.7 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 <typeinfo>
#include <math.h>
#include "sopnamsp.h"
#include "specresp.h"
#include "integ.h"
#include "tvector.h"

/*! 
 * \class SOPHYA::SpectralResponse 
 \ingroup SkyT
 SpectralResponse corresponds to a detector filter's
 response as a function of the frequency.
 The SpectralResponse class is an abstract class.
 The virtual constructor takes as arguments the
 minimum and maximum values of the frequency range on
 which the detector response is defined
 */
//----------------
// Constructor --
//----------------
SpectralResponse::SpectralResponse(double numin, double numax)
{
  _numin = numin;
  _numax = numax;
}


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

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


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

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

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

 
// To change min-max frequency
void
SpectralResponse::setMinMaxFreq(double numin, double numax)
{
  _numin = numin;
  _numax = numax;
}
 

// peakFreq returns the value of the frequency for the
// peak of the spectrum. 
double 
SpectralResponse::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 =transmission(nu);
      if(maxAnswer <= lookForMax) {
	maxAnswer= lookForMax;
	maxNu    = nu;
      }
    }
  return maxNu;    
}


double 
SpectralResponse::peakTransmission()  const
{
  double nuPeak = this->peakFreq();
  return transmission(nuPeak);
}

static SpectralResponse* _mySpecResp = NULL;

static double SpectralResponse_transmission(double nu)
{
   return(_mySpecResp->transmission(nu));   
}


/*! 
  This function performs the integration
  of the transmission function in its frequency range 
  (defined when instanciating an object of on of the
  subclasses)
*/
double 
SpectralResponse::IntegratedSpect()  const
{
  double value = this->IntegratedSpect(_numin, _numax);
  return value;
}

/*! 
  The IntegratedSpect function performs the integration
  of the transmission function in a frequency range <BR> defined by
  numin and numax.
*/
double 
SpectralResponse::IntegratedSpect(double numin, double numax) const 
{
  if(numin < this->minFreq()) numin = this->minFreq();
  if(numax > this->maxFreq()) numax = this->maxFreq();

  
  _mySpecResp = const_cast<SpectralResponse *>(this);
  TrpzInteg I(SpectralResponse_transmission , numin, numax);  
  double val = (double)I;
  _mySpecResp= NULL;
  return(val);
}


static SpectralResponse* _myLogSpecResp = NULL;

static double SpectralResponse_logTransmission(double tau)
{
  double value =  _myLogSpecResp->transmission(pow(10.,tau))*pow(10.,tau);
  return(value);   
}

/*! 
  The logIntegratedSpect function performs the integration
  of the logarithm of the transmission function in a frequency range <BR> defined by
  numin and numax.
*/
double 
SpectralResponse::logIntegratedSpect(double numin, double numax) const 
{
  if(numin <= this->minFreq()) numin = this->minFreq();
  if(numax >= this->maxFreq()) numax = this->maxFreq();
  if(numin == 0) numin = 1.e-99;
  double f1Log = log10(numin);
  double f2Log = log10(numax);
  if(f1Log < -1.e99) f1Log = -1.e99;
  if(f2Log > 1.e99)  f2Log = 1.e99;
  _myLogSpecResp = const_cast<SpectralResponse *>(this); 
  TrpzInteg I(SpectralResponse_logTransmission ,f1Log,f2Log); 
  double val = (double)I;
  _myLogSpecResp= NULL;
  return(val*log(10.));
}

/*! 
  Same than IntegratedSpect(numin,numax) over its
  frequency range defined at the initialisation
*/
double 
SpectralResponse::logIntegratedSpect()  const
{
  double value = this->logIntegratedSpect(_numin, _numax);
  return value;
}


void
SpectralResponse::Print(ostream& os) const 
{

  //   os << "SpectralResponse::Print (" << typeid(*this).name() 
  //<< ") - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
  os << "SpectralResponse::Print - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
  os << "MeanFreq= " << meanFreq() << "  Transmission= " << transmission(meanFreq()) << endl;
  os << "PeakFreq= " << peakFreq() << "  Transmission= " << transmission(peakFreq()) << endl;

}