source: Sophya/trunk/SophyaLib/SkyT/specrespvector.cc@ 2663

//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: specrespvector.cc,v 1.11 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 <math.h>
// #include <typeinfo>

#include "specrespvector.h"
#include "pexceptions.h"
#include "fioarr.h"
//----------------
// Constructor --
//----------------
/*! 
 * \class SOPHYA::SpecRespVec
 \ingroup SkyT
 One may also want to defined the filter of a detector by two vectors:
 one for the frequencies and one for the corresponding transmission values.
 In this case on should use the SpecRespVec class !
*/
SpecRespVec::SpecRespVec()
        : SpectralResponse()
{
}

/*! Constructor: 
  \param nu  vector of the frequencies
  \param fdenu  vector of the corresponding transmission values of the filter
*/
SpecRespVec::SpecRespVec(Vector const & nu, Vector const & fdenu, double numin, double numax)
        : SpectralResponse(numin, numax)
{
  if(nu.NElts() != fdenu.NElts())  
    throw SzMismatchError("SpecRespVec::SpecRespVec() - Non equal vector sizes");
  _vecOfNu = nu;
  _vecOfFDeNu = fdenu;
  _size = nu.NElts();
  if(_vecOfNu.NElts() != _vecOfFDeNu.NElts()) cout << "vectors are not compatible" << exit;
}

SpecRespVec::SpecRespVec(Vector const & nu, Vector const & fdenu)
        : SpectralResponse()
{
  if(nu.NElts() != fdenu.NElts())  
    throw SzMismatchError("SpecRespVec::SpecRespVec() - Non equal vector sizes");
  _vecOfNu = nu;
  _vecOfFDeNu = fdenu;
  _numin = nu(0);
  _numax = nu(nu.NElts()-1);
  _size = nu.NElts();
  if(_vecOfNu.NElts() != _vecOfFDeNu.NElts()) cout << "vectors are not compatible" << exit;
}


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

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



/*! The transmission function extrapolates the transmission values for the 
  frequencies that are not present in the nu vector <BR>
  given at the instanciation of the class.
*/
double 
SpecRespVec::transmission(double nu)  const 
{
  if ( (nu < _numin) || (nu > _numax) ) return(0.);
  double value = 0.;
  int sizeVecOfNu = _vecOfNu.NElts();
  if(nu <=  _vecOfNu(0)) return _vecOfFDeNu(0);
  if(nu >=  _vecOfNu(sizeVecOfNu-1)) return _vecOfFDeNu(sizeVecOfNu-1);

  for (int i=1; i<sizeVecOfNu; i++)
    {
      if(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;
        }
    }
  return value;
}



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

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

}


DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void
ObjFileIO<SpecRespVec>::WriteSelf(POutPersist& s) const
{
  if(dobj == NULL)
    {
      cout << " ObjFileIO<SpecRespVec>::WriteSelf:: dobj= null " << endl;
      return;
    }

  int_4 version, nothing;
  version = 1;
  nothing = 0;   // Reserved for future use 
  s.PutI4(version);
  s.PutI4(nothing);

  s.PutR8(dobj->minFreq());
  s.PutR8(dobj->maxFreq());

  // TVector<T> has Persistence Manager
  s << dobj->getNuVec();
  {
  Vector& xv2 = dobj->getTNuVec();
  cout << xv2 ;
  FIO_TArray<double> vio2(&xv2);
  vio2.Write(s);
  }
}

DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void 
ObjFileIO<SpecRespVec>::ReadSelf(PInPersist& s) 
{
  int_4 version, nothing;
  version = 1;
  nothing = 0;   // Reserved for future use 
  s.GetI4(version);
  s.GetI4(nothing);

  if(dobj == NULL)
    {
    dobj= new SpecRespVec();
    ownobj= true;
    }
  r_8 minf, maxf;
  s.GetR8(minf);
  s.GetR8(maxf);
  dobj->setMinMaxFreq(minf, maxf);
  // TVector<T> has Persistence Manager 
  FIO_TArray<double> vio(&(dobj->getNuVec()));
  vio.Read(s);
  FIO_TArray<double> vio2(&(dobj->getTNuVec()));
  vio2.Read(s);
}


#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template ObjFileIO<SpecRespVec>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class ObjFileIO<SpecRespVec>;
#endif