source: Sophya/trunk/SophyaLib/SkyT/gaussfilt.cc@ 2743

//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: gaussfilt.cc,v 1.9 2004-09-10 09:54:40 cmv Exp $
//
// Description:
//
// History (add to end):
//      Sophie   Oct, 1999  - creation
//
//------------------------------------------------------------------------

//---------------
// C++ Headers --
//---------------
#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <math.h>

#include "gaussfilt.h"

/*! 
 * \class SOPHYA::GaussianFilter
 * \ingroup SkyT
 * Gaussian detector response
 */
GaussianFilter::GaussianFilter()
  : SpectralResponse()
{
  setParams(100, 10., 1.);
}

/*! Constructor: the parameters correspond to the function defined in the 
  <a href="#gausseq"> equation </a> below 
*/
GaussianFilter::GaussianFilter(double nu0, double s, double a, double numin, double numax)
  : SpectralResponse(numin, numax)
{
  setParams(nu0, s, a);
}


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

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

/*! The transmission function is the wel known gaussian:
  <a name="gausseq"> </a>
  \f[
  \hbox{transmission}= A  e^{-{({\nu-\nu_0\over s})^2}};
  \f]
*/

double 
GaussianFilter::transmission(double nu) const
{
  if ((nu < _numin) || (nu > _numax)) return(0.);
  else { 
    double tmp = (nu-_nu0)/_s;
    return(_a * exp(-tmp*tmp));
  }
}

double 
GaussianFilter::peakFreq()  const
{
return(_nu0);
}

double 
GaussianFilter::peakTransmission() const
{
return(_a);
}



void 
GaussianFilter::setParams(double nu0, double s, double a)
{
  if (s < 1.e-19) s = 1.e-19;
  _s = s;
  _nu0 = nu0;
  _a = a;
}



void 
GaussianFilter::Print(ostream& os)  const
{ 
  os << "GaussianFilter::Print - Fmin,Fmax= " << minFreq() << "," << maxFreq() << endl;
  os << " T = A * Exp(-((nu-nu0)/s)^2) : " << " nu0= " << _nu0 << " sig= " 
     << _s << "  A= " << _a << endl;  
  os << "PeakFreq= " << peakFreq() << "  Transmission= " << transmission(peakFreq()) << endl;


}

DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void
ObjFileIO<GaussianFilter>::WriteSelf(POutPersist& s) const
{
  if(dobj == NULL)
    {
      cout << " ObjFileIO<GaussianFilter>::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());
  s.PutR8(dobj->giveNorm());
  s.PutR8(dobj->giveNu0());
  s.PutR8(dobj->giveDNu());
}

DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
void
ObjFileIO<GaussianFilter>::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 GaussianFilter();
      ownobj= true;
    }

  r_8 a, nu0, dnu, numin, numax;
  s.GetR8(numin);
  s.GetR8(numax);
  s.GetR8(a);
  s.GetR8(nu0);
  s.GetR8(dnu);
  dobj->setMinMaxFreq(numin, numax);
  dobj->setParams(nu0, dnu, a);
  //  cout << " Norm - Nu0 - DNu - minFreq - maxFreq " << endl; 
  // cout << _a << "-" << _nu0  << "-" << _s << "-" << _numin << "-" << _numax << endl;
}


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

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