ColdPlasma.h

//-*-C++-*-
/***************************************************************************
 *
 *   Copyright (C) 2006-2010 by Willem van Straten
 *   Licensed under the Academic Free License version 2.1
 *
 ***************************************************************************/
 
// psrchive/More/General/Pulsar/ColdPlasma.h
 
#ifndef __Pulsar_ColdPlasma_h
#define __Pulsar_ColdPlasma_h
 
#include "Pulsar/Transformation.h"
#include "Pulsar/Integration.h"
#include "Pulsar/Archive.h"
#include "Pulsar/Profile.h"
#include "Physical.h"
 
namespace Pulsar {
 
 
  template<class Corrector, class History>
  class ColdPlasma : public Transformation<Integration> {
 
  public:
 
    typedef typename Corrector::Return Type;
 
    ColdPlasma () { name = "ColdPlasma"; }
 
    virtual double get_correction_measure (const Integration*) = 0;
 
    virtual double get_absolute_measure (const Integration*) = 0;
 
    virtual bool get_corrected (const Integration*) = 0;
 
    virtual double get_effective_measure (const Integration*) = 0;
 
    virtual Type get_identity () = 0;
 
    virtual void combine (Type& result, const Type& add) = 0;
 
    virtual void apply (Integration*, unsigned channel) = 0;
 
    virtual void execute (Archive*);
 
    void just_do_it (Archive*);
 
    virtual void revert (Archive*);
 
    void transform (Integration*);
 
 
    void setup (const Integration*);
 
 
    void update (const Integration*);
 
    virtual void set (const Integration* data)
    { setup (data); update (data); }
 
    virtual void set_Profile (const Profile* data)
    { corrector.set_frequency( data->get_centre_frequency () ); }
 
    /* \post All data will be corrected to the reference frequency */
    void execute1 (Integration*);
 
    /* \post All data will be uncorrected to the reference frequency */
    void revert1 (Integration*);
 
    void correct (Integration*, unsigned ichan, unsigned jchan, double freq);
 
    void match (const Integration* reference, Integration* to_be_corrected);
 
    void set_reference_wavelength (double metres);
    double get_reference_wavelength () const;
 
    void set_reference_frequency (double MHz);
    double get_reference_frequency () const;
 
    void set_measure (double measure);
    double get_measure () const;
 
    void set_delta (const Type& d) { delta = d; }
    Type get_delta () const { return delta; }
 
  protected:
 
    friend class Integration;
 
    void range (Integration*, unsigned ichan, unsigned jchan);
 
    Corrector corrector;
 
    Type delta;
 
    std::string name;
 
    std::string val;
 
  private:
    double backup_measure;
  };
 
}
 
 
template<class C, class H>
void Pulsar::ColdPlasma<C,H>::setup (const Integration* data)
{
  set_reference_frequency( data->get_centre_frequency() );
  set_measure( get_correction_measure(data) );
 
  if (Integration::verbose)
    std::cerr << "Pulsar::" + name + "::setup lambda=" 
              << get_reference_wavelength() << " measure=" << get_measure()
              << std::endl;
}
 
template<class C, class H>
void Pulsar::ColdPlasma<C,H>::transform (Integration* data) try
{
  setup (data);
  execute1 (data);
}
catch (Error& error)
{
  throw error += "Pulsar::" + name + "::transform";
}
 
template<class C, class H>
void Pulsar::ColdPlasma<C,H>::execute (Archive* arch)
{
  for (unsigned i=0; i<arch->get_nsubint(); i++)
    execute1( arch->get_Integration(i) );
}
 
template<class C, class H>
void Pulsar::ColdPlasma<C,H>::revert (Archive* arch)
{
  for (unsigned i=0; i<arch->get_nsubint(); i++)
    revert1( arch->get_Integration(i) );
}
 
template<class C, class H>
void Pulsar::ColdPlasma<C,H>::just_do_it (Archive* arch)
{
  delta = get_identity ();
 
  for (unsigned i=0; i<arch->get_nsubint(); i++)
    range (arch->get_Integration(i), 0, arch->get_nchan());
}
 
template<class C, class H>
void Pulsar::ColdPlasma<C,H>::set_reference_frequency (double MHz)
{
  corrector.set_reference_frequency (MHz);
}
 
template<class C, class H>
double Pulsar::ColdPlasma<C,H>::get_reference_frequency () const
{
  return corrector.get_reference_frequency ();
}
 
template<class C, class H>
void Pulsar::ColdPlasma<C,H>::set_reference_wavelength (double metres)
{
  corrector.set_reference_wavelength (metres);
}
 
template<class C, class H>
double Pulsar::ColdPlasma<C,H>::get_reference_wavelength () const
{
  return corrector.get_reference_wavelength ();
}
 
template<class C, class H>
void Pulsar::ColdPlasma<C,H>::set_measure (double measure)
{
  corrector.set_measure (measure);
}
 
template<class C, class H>
double Pulsar::ColdPlasma<C,H>::get_measure () const
{
  return corrector.get_measure ();
}
 
template<class Corrector, class History>
void Pulsar::ColdPlasma<Corrector,History>::update (const Integration* data)
try
{
  backup_measure = get_measure();
 
  double effective_measure = backup_measure;
 
  if (Integration::verbose)
    std::cerr << "Pulsar::" + name + "::update"
                 " backup measure=" << backup_measure << std::endl;
 
  delta = get_identity();
 
  double absolute_measure = get_absolute_measure (data);
  if (absolute_measure)
  {
    if (Integration::verbose)
      std::cerr << "Pulsar::" + name + "::update absolute"
        " measure=" << absolute_measure << std::endl;
 
    // the following corrects centre frequency to infinite frequency
    Corrector absolute;
    absolute.set_measure( absolute_measure );
    absolute.set_reference_wavelength( 0 );
    absolute.set_frequency( data->get_centre_frequency() );
    combine (delta, absolute.evaluate());
 
    // the following corrects channel frequency to centre frequency
    effective_measure += absolute_measure;
  }
 
  const History* corrected = 0;
 
  if ( get_corrected(data) )
    corrected = data->template get<History>();
 
  if (corrected)
  {
    double corrected_measure = corrected->get_measure();
    double lambda = corrected->get_reference_wavelength();
 
    if (Integration::verbose)
      std::cerr << "Pulsar::" + name + "::update corrected"
        " measure=" << corrected_measure << " lambda=" << lambda 
                << std::endl;
 
    // calculate the correction due to the new centre frequency, if any
    corrector.set_wavelength( lambda );
    combine (delta, corrector.evaluate());
 
    effective_measure -= corrected_measure;
  }
 
  if (Integration::verbose)
    std::cerr << "Pulsar::" + name + "::update effective_measure="
              << effective_measure << std::endl;
 
  // set the effective correction measure
  set_measure( effective_measure );
}
catch (Error& error)
{
  throw error += "Pulsar::"+name+"::update";
}
 
template<class C, class History>
void Pulsar::ColdPlasma<C,History>::execute1 (Integration* data) try
{
  update (data);
 
  if (Integration::verbose)
    std::cerr << "Pulsar::"+name+"::execute1"
      " effective "+val+"=" << get_measure() <<
      " reference wavelength=" << get_reference_wavelength() << std::endl;
 
  range (data, 0, data->get_nchan());
 
  // restore the original correction measure
  set_measure( backup_measure );
 
  History* corrected = new History;
  corrected->set_measure( backup_measure );
  corrected->set_reference_wavelength( get_reference_wavelength() );
  data->add_extension( corrected );
 
}
catch (Error& error)
{
  throw error += "Pulsar::"+name+"::execute1";
}
 
template<class C, class History>
void Pulsar::ColdPlasma<C,History>::revert1 (Integration* data) try
{
  if ( !get_corrected(data) )
    return;
 
  History* corrected = data->template get<History>();
 
  if (!corrected)
    throw Error (InvalidState, "Pulsar::" + name + "::revert1",
                 "no correction history");
 
  set_measure( -corrected->get_measure() );
  set_reference_wavelength( corrected->get_reference_wavelength() );
  set_delta( get_identity() );
 
  if (Integration::verbose)
    std::cerr << "Pulsar::"+name+"::revert1"
      " effective "+val+"=" << get_measure() <<
      " reference wavelength=" << get_reference_wavelength() << std::endl;
 
  range (data, 0, data->get_nchan());
 
  // this should remove the history
  delete corrected;
}
catch (Error& error)
{
  throw error += "Pulsar::"+name+"::revert1";
}
 
template<class C, class History>
void Pulsar::ColdPlasma<C,History>::match (const Integration* reference,
                                           Integration* to_correct)
{
  const History* corrected = reference->template get<History>();
  if (!corrected)
    throw Error (InvalidState, "Pulsar::" + name + "::match",
                 "reference has no correction history");
 
  set_measure( corrected->get_measure() );
  set_reference_wavelength( corrected->get_reference_wavelength() );
 
  execute1( to_correct );
}
 
template<class C, class H>
void Pulsar::ColdPlasma<C,H>::range (Integration* data,
                                     unsigned ichan, unsigned kchan) try
{
 
  if (Integration::verbose)
    std::cerr << "Pulsar::"+name+"::range "+val+"=" << get_measure()
         << " lambda_0=" << get_reference_wavelength() << " m" 
         << " delta=" << delta << std::endl;
 
  if (get_measure() == 0 && delta == get_identity())
    return;
 
  if (ichan >= data->get_nchan())
    throw Error (InvalidRange, "Pulsar::"+name+"::range",
                 "start chan=%d >= nchan=%d", ichan, data->get_nchan());
 
  if (kchan > data->get_nchan())
    throw Error (InvalidRange, "Pulsar::"+name+"::range",
                 "end chan=%d > nchan=%d", kchan, data->get_nchan());
 
  for (unsigned jchan=ichan; jchan < kchan; jchan++)
  {
    corrector.set_frequency( data->get_centre_frequency (jchan) );
    apply (data, jchan);
  }
 
}
catch (Error& error) {
  throw error += "Pulsar::"+name+"::range";
}
 
 
template<class C, class H>
void Pulsar::ColdPlasma<C,H>::correct (Integration* subint,
                                       unsigned ichan, unsigned kchan, 
                                       double reference_frequency) 
try
{
  double effective_measure = get_effective_measure (subint);
 
  if (Integration::verbose)
    std::cerr << "Pulsar::" + name + "::correct"
                 " reference freq=" << reference_frequency << 
                 " effective measure=" << effective_measure << std::endl;
 
  set_measure( effective_measure );
  set_reference_frequency( reference_frequency );
  set_delta( get_identity() );
 
  range (subint, ichan, kchan);
}
catch (Error& error)
{
  throw error += "Pulsar::"+name+"::correct";
}
 
 
#endif