source: Sophya/trunk/SophyaLib/SkyT/radspecvector.cc@ 601

//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: radspecvector.cc,v 1.1.1.1 1999-11-19 16:34:32 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 <math.h>
// #include <typeinfo>

#include "radspecvector.h"
#include "pexceptions.h"

//----------------
// Constructor --
//----------------
RadSpectraVec::RadSpectraVec(Vector const & nu, Vector const & fdenu, double numin, double numax)
        : RadSpectra(numin, numax)
{
  if(nu.NElts() != fdenu.NElts())  
    throw SzMismatchError("RadSpectraVec::RadSpectraVec() - Non equal vector sizes");
  _vecOfNu = nu;
  _vecOfFDeNu = fdenu;
  _numin = nu(0);
  _numax = nu(nu.NElts()-1);
  if(_vecOfNu.NElts() != _vecOfFDeNu.NElts()) cout << "vectors are not compatible" << exit;
}


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

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




double 
RadSpectraVec::flux(double nu) const
{
  if ( (nu < _numin) || (nu > _numax) ) return(0.);
  double value = -99;
  int sizeVecOfNu = _vecOfNu.NElts();
  double minVal = _vecOfNu(0);
  if(nu <= minVal) return _vecOfFDeNu(0);
  if(nu >=  _vecOfNu(sizeVecOfNu-1)) return _vecOfFDeNu(sizeVecOfNu-1);
  
  for (int i=1; i<sizeVecOfNu; i++)
    {
      if(nu>= minVal && nu < _vecOfNu(i))
        {
          double up = _vecOfFDeNu(i) ;
          double down = _vecOfFDeNu(i-1);
          double xmin = _vecOfNu(i-1);
          double xmax = _vecOfNu(i);
          double a = ((up-down)/(xmax-xmin));
          value = a*nu+(up-a*xmax);
          return value;
        }
      else 
        {
          minVal = _vecOfNu(i);
        }       
    }
  return value;
}


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

}