source: Sophya/trunk/SophyaLib/SkyT/specresp.cc@ 610

//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: specresp.cc,v 1.2 1999-11-20 21:00:52 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 <typeinfo>
#include <math.h>
#include "specresp.h"
#include "integ.h"

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

 

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


double 
SpectralResponse::IntegratedSpect()  const
{
  double value = this->IntegratedSpect(_numin, _numax);
  return value;
}

double 
SpectralResponse::IntegratedSpect(double numin, double numax) const 
{
  _mySpecResp = const_cast<SpectralResponse *>(this);
  if(numin <= _numin) numin = _numin;
  if(numax >= _numax) numax = _numax;
  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);   
}

double 
SpectralResponse::logIntegratedSpect(double numin, double numax) const 
{
  if(numin <= _numin) numin = _numin;
  if(numax >= _numax) numax = _numax;
  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.));
}

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;

}